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1 3.  ABSTRACT  ( Maximum  200  words) 

An  abstract,  algebraic  bilevel  version  of  conventional  multigrid  methods  has  been  developed  that  formally  unifies  and  extends  the  reduced  basis, 
substructuring,  smoothing/homogenization,  and  frequency  window  methods  of  model  reduction.  Within  the  context  of  this  formalism,  a  given  model 
reduction  method  synthesizes  a  reduced  model  that  plays  the  role  of  a  “coarse  grid”  model.  Another  original-model  approximation,  conjugate  to  the 
model  reduction  method,  plays  the  role  of  a  “fine  grid”  model.  Conventional  multigrid  methods  can  be  thought  of  as  an  extension  of  the  coarse  grid 
model  beyond  the  original  scope  of  its  accuracy  through  the  alternating  use  of “fme-grid-error-smoothing”  and  “coarse -grid-correction”  steps.  Analo¬ 
gously  fashioned,  abstract  versions  of  these  steps  extend  the  particular  model  reduction  method  used  beyond  the  original  scope  of  its  accuracy.  In 
addition  to  the  development  of  the  mathematical  formalism,  a  generic  continuation  is  proposed  and  developed  as  the  conjugate  approximation  for  use 
in  the  abstract  version  of  the  fme-grid-error-smoothing  step.  For  generality,  the  nature  of  the  parameter-embedding  of  the  continuation  is  deliberately 
left  unspecified;  a  particular  parameter-embedding  is  chosen  by  the  analyst  to  complement  the  particular  model  and  model  reduction  method  (or  class 
thereof)  in  a  given  case.  As  an  example,  a  particular  submethod  of  this  formalism  is  constructed  by  specializing  the  parameter-embedding  to  a 
temporal,  biscale  perturbational  parameter-embedding  for  the  case  of  a  generic  set  of  coupled  ordinary  differential  equations  with  constant  coeffi¬ 
cients.  The  initial  iteration  of  this  submethod  is  a  general,  combined  transient/frequency-window  methodology  for  linear  finite -element  method 
applications.  It  is  shown  to  encompass  both  the  force  derivative  and  Lanczos/Ritz  vector  methods  for  the  limiting  case  of  a  zero  frequency,  single¬ 
point  window  for  the  “fast”  time  scale.  To  put  this  in  perspective,  the  force  derivative  and  Lanczos/Ritz  vector  methods  are  two  leading  approaches  for 
extending/replacing  modal  reduced  basis  methods.  Current  frequency  window  methods,  in  turn,  are  very  efficient  for  extensive  time-harmonic 
reanalysis  of  the  system. 
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A  GENERIC  BILEVEL  FORMALISM  FOR  UNIFYING 
AND  EXTENDING  MODEL  REDUCTION  METHODS 


1.  INTRODUCTION 

In  its  broadest  sense,  model  reduction  methods  synthesize,  from  one  model,  a  new  “reduced”  model 
whose  response-predictions  coincide  with  that  of  the  original  model  for  a  specified  subset  of  the  original 
scope  of  stimuli.  Within  this  restricted  scope  of  admissible  stimuli  driving  the  response,  the  synthesized 
model  can  be  used  in  place  of  the  original  model.  The  benefit  usually  corresponding  to  this  is  a  considerable 
reduction  in  overall  computational  effort.  This  is  satisfactory  if  the  modeler/analyst  is  interested  only  in  the 
model’s  response  to  this  subset  of  admissible  stimuli.  The  problem  remains,  however,  that  outside  of  this  re¬ 
stricted  scope  of  admissible  stimuli,  the  synthesized  model’s  predictions  typically  compare  poorly  with  that 
of  the  original  model.  A  remedy  to  this  dramatic  dropoff  in  accuracy  is  to  somehow  approximate  this  “resid¬ 
ual”  part  of  the  model  in  a  manner  that  is  consistent  with  the  reduced-model  synthesized  and,  at  the  same 
time,  is  complementary  to  the  model  reduction  method.  This  results  in  a  hybrid  approximation  to  the  origi¬ 
nal  model.  “Consistent”  here  means  that  the  scope  over  which  the  reduced  model  is  faithful  to  the  original 
model  is  fully  preserved  when  it  is  combined  with  the  other  approximation,  conjugate  to  the  model  reduction 
method,  to  form  the  hybrid.  “Complementary”  here  means  that  the  conjugate  approximation  be  an  adequate 
one  over  at  least  the  set  of  stimuli  not  admissible  to  the  corresponding  model  reduction  method.  A  good 
role  model  for  the  combination  of  complementary  approximations  into  a  self-consistent,  efficient  hybrid  is 
that  of  multigrid  methods  [1,2].  One  approximation  is  effective  at  eliminating  oscillatory  error  components 
(on  a  fine  grid)  and  the  other  is  effective  at  eliminating  smooth  error  components  (on  a  coarse  grid).  In  this 
paper,  an  abstract,  algebraic  bilevel  version  of  this  approach  is  developed  that  provides  a  general  means  of 
consistently  combining  model  reduction  methods  with  other  approximations.  The  generic  formalism  en¬ 
compasses  at  least  three  broad  classes  of  model  reduction  methods:  smoothing/homogenization,  reduced 
basis,  and  substructuring  methods.  A  particular  model  reduction  method  operates  at  one  “level”  (the  coarse 
grid  analogue)  and  some  conjugate  approximation  (of  the  original  model  response)  operates  at  the  other 
“level”  (the  fine  grid  analogue).  A  generic  continuation  approach  is  proposed  and  developed  in  this  paper 
as  a  class  of  conjugate  approximations.  This  choice  is  seen  to  unify  and  generalize  several  successful  ap¬ 
proaches  to  extending  the  utility  of  reduced  basis  methods. 

One  important  class  of  model  reduction  methods  applies  to  the  modeling  of  coupled-multiscale  phenom¬ 
ena  for  which  the  differences  in  scale  are  significant.  A  prominant  example  of  this  would  be  a  structural- 
component  model  of  a  material  that  possess  an  intricate  spatial  heterogeneity  of  a  length  scale  small  with 
respect  to  structural  dimensions  but  large  with  respect  to  atomic  dimensions.  Homogenization  and  smooth¬ 
ing  methods  of  model  reduction  were  developed  for  these  and  similar  such  cases.  Both  methods  refer  to 
the  process  of  mathematically  synthesizing  a  macroscale  (usually  effective  single-phase  constitutive)  model 
from  a  given  mesoscale/microscale  model  such  that  the  predictions  of  each  coincide  on  the  macroscale. 
Methods  designed  for  periodic  media  are  usually  referred  to  as  homogenization  methods  in  the  literature. 

Manuscript  approved  September  30,  1999. 
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The  best  known  and  most  widely  used  is  the  asymptotic  approach  [3,  4,  5,  6, 7],  which  is  based  on  special¬ 
izing  the  derivative-expansion  version  Ref  [8,  pp.  230-232]  of  the  multiple  scales  perturbation  method  to 
a  biscale  expansion  in  the  spatial  variables.  Methods  of  smoothing  (see  Refs.  [9]  and  [10,  p.  16]),  which 
are  designed  for  statistically  homogeneous  random  media,  have  several  variants,  including  “self-consistent” 
methods  [10,  11]  and  methods  that  bound  [10,  12]  the  constitutive  parameter  values  through  the  use  of 
variational  principles.  Fishman  and  McCoy  [13]  showed  that  both  homogenization  and  smoothing  meth¬ 
ods  can  be  united  under  a  single  projection  operator  framework.  Similar  projection  techniques  have  been 
used  previously  to  deduce  the  influence  of  one  scale  upon  another  in  other  areas,  such  as  neutron  transport 
Refs.  [14,  pp.  37,  193-194,  and  236-250]  and  nonequilibrium  statistical  mechanics  [15].  An  example  of  a 
type  of  heterogeneity  that  does  not  fit  either  of  the  usually-assumed  periodic  or  random  material  models  is 
found  in  Steinberg  and  McCoy  [16].  They  consider  the  case  of  fluid-loaded  stmctures  for  which  the  hetero¬ 
geneous  regions,  each  of  limited  spatial  extent  on  the  macroscale,  are  irregularly  distributed.  The  enabling 
assumption  of  the  projection  technique  [13, 16]  underlying  both  homogenization  and  smoothing  methods  is 
that  the  fluctuation-scale-response  is  excited  solely  by  its  coupling  to  the  macroscale  response;  any  external 
stimulation  of  the  system  is  taken  to  occur  at  the  macroscale  only. 


Multiscale  systems  are  also  represented  by  models  with  a  finite  number  of  degrees-of-freedom  (DOF). 
If  properly  chosen,  a  given  subset  of  the  model’s  DOF  can  represent  a  corresponding  scale  of  the  model, 
with  fewer  DOF  usually  required  the  larger  the  scale.  Conventional  multigrid  methods  can  be  used  to  model 
these  scales  and  their  coupling  for  models  consisting  of  sets  of  algebraic  equations.  Another  important  class 
of  finite-degree-of -freedom  mathematical  models  consists  of  coupled  sets  of  ordinary  differential  equations, 
a  common  source  of  such  models  being  spatial  discretization  of  a  physical  system  via  finite-element  meth¬ 
ods  (FEMs)  or  finite-difference  methods.  The  number  of  DOF  of  such  models  is  often  large,  making  their 
response-prediction  computationally  expensive.  Two  overlapping  classes  of  DOF  reduction  methods  appro¬ 
priate  to  such  systems  are  reduced  basis  methods,  which  use  a  Rayleigh-Ritz  approximation  with  respect  to 
a  reduced  set  of  generalized  coordinates,  and  “substructuring”  methods,  for  which  the  DOF  of  the  synthe¬ 
sized  reduced-model  consists  of  a  chosen  subset  of  the  DOF  of  the  original  model.  The  admissible  stimuli 
for  the  reduced  basis  method  consists  of  the  span  of  a  small  number  of  basis  vectors  (the  “reduced  basis”), 
the  number  being  small  compared  with  the  total  number  of  DOF  if  the  method  is  to  be  efficient.  In  turn, 
the  usual  enabling  assumption  for  the  substructuring  methods  is  that  the  eliminated  DOF,  that  is,  those  not 
retained  in  the  reduced-model  (and  typically  composing  the  “substructures”),  cannot  be  externally  stimu¬ 
lated  (loaded).  A  partial  summary  of  substructuring  methods  is  given  by  Abdelhamid  [17].  The  first  and 
still  very  popular  such  method  is  that  due  to  Guyan  [18].  Flippen  [19]  derived  a  number  of  such  methods 
from  a  general  time-derivative  series  solution  [20].  The  Modal  Reduction  method  [21,  22,  23],  popular  for 
the  synthesis  of  Test-Analysis-Models  (TAMs)  via  elimination  of  (almost)  all  but  the  test-sensor-associated 
DOF,  is  an  adaptation  of  modal  reduced  basis  methods  to  substructuring. 


For  reduced  basis  methods,  the  reduced  basis  sets  consist  mainly  of  either  Lanczos/Ritz,  modal,  or 
solution-path-derivative  vectors.  As  pointed  out  by  Nour-Omid  and  Clough  (Ref.  [24, 'p.  566]),  Ritz  vector 
methods  [25],  as  usually  implemented,  are  essentially  Lanczos  methods.  Lanczos  methods  are  applicable  to 
models  constrained  to  loads  contained  within  the  span  of  a  small  number  of  fixed  vectors.  These  fixed  load- 
basis  vectors  are  used,  in  turn,  to  construct  a  reduced  basis  spanning  the  Krylov  subspace  associated  with 
the  method.  The  Lanczos  method  is  usually  discussed  in  terms  of  a  second-order  formulation,  but  it  also 
has  a  first-order  formulation  [26]  for  nonproportional  damping,  as  well  as  block  [27],  modal-hybrid  [28], 
and  other  variants.  Noting  the  analogy  of  their  method  to  those  based  on  Krylov  subspaces  (in  particular, 
those  of  Amoldi  and  Lanczos),  Haggblad  and  Eriksson  [29]  recognized  that  efficient,  “generalized  Krylov” 
subspace  descriptions  can  be  generated  as  recursive  relations  from  series  solutions  of  the  governing  equa- 
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tions.  They  take  a  physical  approach  and  derive  advanced  methods  of  this  type  from  series  in  which  “each 
component  is  successively  computed  by  balancing  the  inertia  forces  from  the  previous  component  in  the  se¬ 
ries.”  This  prescription  for  generating  their  series  is  model-specific  to  mechanics.  In  particular,  they  do  not 
explicitly  identify  an  underlying  generic  perturbation  or  continuation  method  associated  with  any  of  their 
series.  The  application  of  modal  reduced  basis  methods  for  linear  problems  is  commonplace.  Its  applica¬ 
tion  has  been  extended  [30]  to  nonlinear  structural  dynamics  problems  as  well.  Continuation-based  reduced 
basis  methods,  as  summarized  mathematically  by  Rheinboldt  (Ref.  [31,  pp.  80-86])  for  example,  seem  to 
have  a  growing  history  of  practical  success  with  respect  to  some  nonlinear  problems  (see  Noor  [32],  Noor 
and  Whitworth  (Ref.  [33,  p.  915]),  Noor  (Ref.  [34,  p.  958]),  and  Noor  and  Peters  [35],  for  example).  The 
solution  and  its  path  derivatives  at  a  point  on  the  continuation  path  are  used  as  a  Rayleigh-Ritz  reduced  basis. 

Kammer  [21]  recognized  the  need  for  and  constructed  a  hybrid  extension  to  modal  reduction  for  build¬ 
ing  robust  TAMs.  The  force  derivative  method  is  another  example  of  a  model  reduction  hybrid,  in  this  case 
for  extending  modal  reduced  basis  approximations.  With  linear  FEM  models  of  mechanical  systems  as  a 
benchmark,  comparisons  [36,  37,  38]  between  several  reduced  basis  methods  based  on  computed  responses 
of  the  resulting  reduced  dynamical  models,  seem  to  favor  the  force-derivative  [37,  39,  38]  and  the  (closely 
related)  Lanczos  [24,  27,  26,  40]  methods.  This  is  evidence  supporting  the  utility  of  hybrid  extensions  to 
model  reduction  methods.  The  force-derivative  method  encompasses  [39,  38]  the  older  mode-displacement 
and  mode-acceleration  methods  [41,  42]  as  zeroth  and  first-order  submethods,  respectively.  It  is  essen¬ 
tially  a  modal  basis  method  that  is  systematically  corrected  by  terms  containing  residual  mode  information. 
(The  chosen  modal  basis  vectors  are  the  target  modes  whose  span  contains  the  reduced  model’s  response.) 
Its  importance  is  underscored  by  the  fact  that  it  provides  both  a  foundation  and  extension  to  the  popular 
Craig-Bampton  [43]  version  of  component  mode  synthesis  [44],  as  shown  by  Suarez  and  Singh  [45].  The 
conventional  version  of  the  force  derivative  method  is  a  convolution  integral  formulation  [39,  38],  from 
which  higher  order  corrections  are  derived  via  repeated  integration  by  parts  of  the  integral.  This  integration- 
by-parts  development  of  the  series  is  apparently  serendipitous;  a  direct  systematic  development  in  terms 
of  a  generic  perturbation  methodology  is  neither  mentioned  nor  used  as  an  alternative.  As  shown  later  in 
this  paper,  the  force  derivative  and  Lanczos  methods  are  closely  related.  For  the  case  of  generic  undamped 
mechanical  FEM  models,  if  the  force  derivative  method  is  constrained  to  loads  within  the  span  of  a  small 
number  of  fixed  vectors  (a  Lanczos  requirement),  each  correction  term  of  the  method  then  reduces  to  a 
scalar-multiple  of  a  basis  vector  for  the  (Lanczos-associated)  Krylov  subspace. 

Frequency  window  methods  form  another  class  of  hybrid  methods,  their  original  intention  essentially 
being  an  extension  to  modal  methods.  They  are  very  efficient  for  extensive  time-harmonic  re-analysis.  In 
structural  acoustics,  for  example,  fine  frequency  sweeps  may  be  needed  to  build  transient  responses  (Ref  [46, 
p.  251])  using  fast  Fourier  transforms.  Igusa  and  Achenbach  et  al.  [47,  48,  46]  have  developed  frequency 
window  methods  for  which  substructure  attachments  are  coupled  to  a  main  structure  by  Lagrange  multipli¬ 
ers.  Its  efficiency  derives  from  its  use  of  two  complementary  approximations,  a  frequency -response  repre¬ 
sentation  of  the  resonances  by  simple  analytical  forms  in  conjunction  with  frequency-interpolation  over  the 
nonresonance  part  of  the  response  within  a  “window.”  In  Ref.  [46],  for  example,  the  modes  of  a  fluid-loaded 
shell  are  used  to  analyze  the  response  of  the  same  shell  with  internal  substructures.  The  eigensolution  is 
“subtracted  out”  of  the  fluid-shell  response,  leaving  an  interpolation  over  the  “smooth,”  nonresonance  part 
of  the  fluid-shell  response  within  a  window.  An  analytical  expression  is  used  to  represent  the  resonant  part 
of  the  fluid-shell  response  in  terms  of  the  the  fluid-shell  eigenvalues  and  eigenvectors,  which  are  obtained  by 
independent,  external  means.  Flippen  [49]  developed  a  frequency  window  version  of  substructuring  meth¬ 
ods  for  degree-of-freedom  reduction.  Ingel  et  al.  [50]  extended  this  into  a  finite-element  environment. 
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The  formalism  of  this  paper  attempts  to  encompass  these  currently  accepted,  in-use  hybrid  methods  as 
special  limiting  cases  while  mathematically  extending  their  foundation  to  a  more  general  model  reduction 
setting.  In  this  context,  specific  computational  implementations  and  applications  are  beyond  the  scope  of 
this  report.  Feedback  from  such  computational  experience  is  ultimately  necessary  for  refining  an  algoritm 
and  increasing  its  efficiency.  Therefore,  the  practical  envelope-of-utility  of  any  particular  algorithm  arising 
from  this  formalism  cannot  be  fully  evaluated  here. 

2.  MATHEMATICAL  FORMALISM 

Let  the  mathematical  model  to  be  considered  in  this  report  be  generically  represented  by 

Lu  =  /,  (2.1) 

where  L  is  the  system  operator,  u  is  the  system  response,  and  /  is  the  system  stimulus  (loads,  sources, 
etc.)  driving  the  system.  The  formal,  generic  nature  of  the  model  description  (2.1)  allows  for  the  derivation 
of  results  of  a  general  nature  and  utility.  It  should  be  noted  that  (2.1)  is  not  necessarily  in  a  state-space 
formulation.  Although  some  of  the  results  that  follow  are  not  new,  they  are  included  for  completeness.  In 
addition,  the  development  that  follows  carries  with  it  the  implicit  caveat  that  the  ranges  and  domains  of 
the  appropriate  operators  are  such  that  the  various  operator  compositions  indicated  are  well-defined.  As  an 
algebraic  notation,  juxtapositioning  of  operators  in  this  report  denotes  their  compositioning  as  mappings. 
The  additional  notation  7v(  A)  =  {range  of  A},  V(  A)  =  {domain  of  A},  and  Af  (A)  =  {nullspace  of  A}  for 
generic  operator  A  will  also  be  used. 

2.1  Generalized  Inverse  Theory  for  Model  Reduction 

The  typical  goal  of  any  particular  model  reduction  method  is  to  provide  an  accurate  approximation  to  L~l 
over  the  “relevant”  subsets  of  possible  u9 s  and  /’ s  of  (2.1),  that  is,  over  a  restricted  domain  and  range  for  L. 
In  this  report,  this  approximation  to  L~l  is  made  generic  to  all  such  model  reduction  methods  and  system 
models  (2.1)  through  the  concept  of  an  “outer  generalized  inverse”  of  L . 

Definition  1  For  a  given  operator  L,  an  outer  generalized  inverse  of  L,  denoted  by  L1 ,  satisfies 

L1 LL 1  =  L1.  (2.2) 

Similarly,  an  “inner  generalized  inverse”  of  L,  denoted  by  Ln ,  satisfies 


LLnL  =  L.  (2.3) 

The  ordinary  inverse  L~1,  when  it  exists,  is  both  an  inner  and  outer  generalized  inverse.  This  report  is 
exclusively  concerned  with  outer  generalized  inverses.  The  terminology  is  borrowed  from  the  theory  of  ma¬ 
trices  (Ref.  [51,  pp.  428^432]),  for  which  a  matrix  generalized  inverse  satifying  both  (2.2)  and  (2.3)  always 
exists  for  any  given  matrix  L.  In  a  more  general  operator  setting  for  which  L  is  not  necessarily  a  matrix, 
one  might  satisfy  one  of  either  (2.2)  or  (2.3)  without  satisfying  the  other.  The  distinction  between  which  of 
(2.2)  and  (2.3)  is  satisfied  is  then  necessary.  The  terms  “inner”  and  “outer”  are  used  to  make  this  distinction 
for  lack  of  a  known  precedence  regarding  terminology.  As  in  the  matrix  case,  generalized  inverses  in  this 
more  general  setting  are  not,  in  general,  unique.  The  following  theorem  also  carries  over  from  the  matrix 
case  to  the  general  operator  setting. 
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Theorem  1  The  operators  LL1  and  L1  L  are  each  idempotent.  If  L  is  bijective,  so  that  L  1  exists,  then 

L-1  =  Lr  +  (I-L1L)L~1  (2.4) 

L"1  =  L1  +  L~l(I  -  LL1)  (2.5) 

are  identities  expressing  L  ~ 1  in  terms  of  L 1  and  a  “residual”  term  containing  one  of  the  idempotent  operators 

(/ -  LiL)ot{I  -  LL1). 

The  proof  that  LL1  and  L1  L  are  each  idempotent  follows  directly  from  (2.2).  (An  idempotent  .4  satisfies 
.4  2  =  .4.4  =  .4.)  The  proof  of  (2.4)  and  (2.5)  follows  from  the  decompositons  of  the  identity  given  by 

I  =  L1  L  +  (/  -  L1  L)  (2.6) 

/  =  LL1  +  [I  -  LL1),  (2.7) 

respectively.  In  (2.6)  each  element  of  V  ( L )  is  uniquely  decomposed  into  the  sum  of  a  component  in  7v  ( Z/i ) 
and  a  component  in  TZ(I  -  L1  L).  Similarly,  in  (2.7)  each  element  of  7 Z{L)  is  uniquely  decomposed  into 
the  sum  of  a  component  in  LULL1)  and  a  component  in  7v(/  -  LL1).  If  one  uses  L 1  as  an  approximation 
to  I/-1  in  (2.4)  or  (2.5)  by  neglecting  the  residual  term  (for  now),  one  is  effectively  computing  a  reduced  set 
of  responses  LI {LrL)  for  (2.1)  for  a  reduced  set  of  stimuli  LKLL1).  This  can  be  more  easily  seen  by  use  of 
the  identities 

L1  =  L-^LL1)  (2.8) 

L1  =  (LIL)L~l.  (2.9) 

The  first  identity  shows  that  filtering  the  stimuli  in  7 Z{L)  by  LL1  and  then  solving  (2.1)  is  equivalent  to 
using  L1  directly.  Similarly,  the  second  identity  shows  that  filtering  the  response  (obtained  from  solving 
(2.1))  by  L1  L  is  also  equivalent  to  using  L1  directly. 

To  make  use  of  outer  generalized  inverses  as  model  reduction  approximations,  one  must  be  able  to  con¬ 
struct  them  so  that  either  TZ{LLt),  7 Z(LrL),  or  both,  are  adjusted  to  correspond  to  the  set  of  relevant  stim¬ 
uli,  responses,  or  both,  respectively.  The  following  theorem  gives  conditions  for  predetermining  7 Z{LL!), 
7v(L/L),  or  both  for  a  certain  class  of  outer  generalized  inverses. 


Theorem  2  Given  a  system  operator  L  for  system  model  (2.1),  let  the  relevant  stimuli  and  response  subsets 
be  given  by  LZ{Pr)  C  LZ{L)  and  lZ{Pd)  C  T>(L),  respectively,  for  the  linear  idempotent  operators  Pr  and 
Pj.  Define  the  effective  version  of  L  as 


Lcjj  =  P-LPi. 

(2.10) 

If  I7"  is  defined  by 

L1*  =  PjLrj/Pr. 

(2.11) 

then  L  /x  is  an  outer  generalized  inverse  of  L  and 

X 

1 

Ul 

cc 

(2.12) 

For  a  given  /  in  (2.1),  define  a  model  reduction  approximation  to  the  response  as  u 
in  (2.11)  satisfies  either 

=  L1"  f  for  which  Ltjj 1 

*-* 

II 

a? 

(2.13) 

or 

(2.14) 

or  both. 
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•  If  Lej j1  satisfies  (2.13)  then 


TZ{Lt*L) 

Lr*LPd 

men) 

•  If  Lejf1  satisfies  (2.14)  then 

7Z(LLr) 

PrLL1* 

N{PrL) 


n(Pd) 

A'(  /  -  L1* L)  (2.15) 

Pd  (2.16) 

Ar(Pd).  (2.17) 

K(Pr) 

A'(7  -  LL1*)  (2.18) 

Pr  (2.19) 

7Z{I  -  L1*  L).  (2.20) 


•  If  satisfies  both  (2.13)  and  (2.14),  then  L1*  is  unique. 

•  If  Pr  and  Pd  satisfy 

PrL{I-Pd)  =0 

and  Lzfj1  satisfies  (2.13),  then 

LIXL  =  Pd 


and  Lej  j1  satisfies  (2.14)  as  well. 

•  If  Pr  and  P,j  satisfy 

(I-Pr)LPd  =  0 

and  Lef  j1  satisfies  (2.14),  then 

LL1  =Pr. 


If,  in  addition,  Ar{L)  =  {0}  then  Lejj!  satisfies  (2.13)  as  well. 


(2.21) 

(2.22) 


(2.23) 

(2.24) 


The  subscripts  “r”  and  “d”  on  Pr  and  Pd  denote  range  and  domain  (of  L),  respectively.  Appendix  A  con¬ 
tains  the  proof  of  Theorem  2.  Theorem  2  reduces  the  model  reduction  process  to  finding  a  L,  j  j 1  satisfying 
either  (2.13)  to  obtain  TZ(LrL)  =  7 Z[Pd)  of  (2.15),  (2.14)  to  obtain  lZ(LLr*)  =  7 Z(Pr)  of  (2.18),  or  both. 
The  L.  jf  approximation  of  L  given  by  (2.10)  is  essentially  a  generalized  Petrov  (or  Galerkin)  method  for 
arbitrary  (not  necessarily  of  finite  dimension)  subspaces. 

A  primary  advantage  of  satisfying  either  of  the  constraints  (2.21)  or  (2.23)  is  that  L1*  f  may  then  be  an 
exact  solution  to  (2.1).  In  the  case  for  which  (2.22)  is  valid,  if  u  =  L!~f,  then  Pdu  =  Lr'Lu  =  L^LL1*  f  = 
L1  f  =  u  and  L1  (Lit  —  f)  =  Pdu  —  L1*  f  =  Pdu  —  u  =  0,  so  that  (Lu  —  f)  €  A '(L1*).  In  the  case  for 
which  (2.24)  is  valid,  /  =  Pr  f  =  LL1  f  =  Lu  shows  that  u  =  L1*  f  is  an  exact  solution  to  (2.1)  for  all  / 
satisfying  Pr  f  =  f. 

The  constraint  (2.23)  on  Pr  and  Pd  for  model  reduction  is  not  new.  In  fact,  PrLPd  =  LPd  corresponds 
exactly  with  the  constraint  (6)  of  (Ref.  [52,  p.  126])  when  making  the  associations  Pr  P  and  Pd  QP, 
where  P  and  9.P  are  the  notations  of  (Ref.  [52]).  The  constraint  (3)  of  (Ref.  [52,  p.  125])  on  stimuli  / 
justifies  the  association  Pr  -4  P.  The  PrLPd  =  LPd  constraint,  in  conjunction  with  (2.14),  is  seen  here  to 
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arise  naturally  as  a  means  of  simultaneously  obtaining  both  (2.24)  and  1Z(LJ  L)  =  l\  ( If )  of  (2.15)  in  the 
case  of  bijective  L.  A  similar  advantage  is  acmed  by  adherence  to  (2.13). 

The  following  corollary  of  Theorem  2  provides  alternatives  to  (2.4)  and  (2.5)  as  decompositions  of  T_1 
into  Lr'  and  a  residual  term. 

Corollary  1  If,  in  Theorem  2,  L  is  bijective  then  (2.13)  implies 

L~l  =  Lr  +  (I-  LrL)(I-  Pd)L-]  (2.25) 

and  (2.14)  implies 

L-1  =  L1*  +  L~1(I  -  Pr){I  -  LL1*).  (2.26) 


Proof:  The  relation  (2.13)  implies  (2.16),  which  implies  (/  -  L1* L)Pd  =  0,  which  in  turn  implies  ( I  - 
L1* L){I  -  Pd)  =  (I  -  LrL),  and  (2.25)  follows  from  (2.4).  The  relation  (2.14)  implies  (2.19),  which 
implies  Pr{I  —  LL1*)  =  0,  which  in  turn  implies  (I  -  Pr)(7  —  LL1  )  =  (I  —  LL 1  ),  and  (2.26)  follows 
from  (2.5). 

2.2  A  Bilevel  Formalism 

Taking  L1  — >  L 1  in  (2.4)  and  operating  upon  /,  with  L~]  j  ~  Uj+i  on  the  left-hand  side  and  l~lf  ~  uj 
on  the  right-hand  side,  suggests  the  iteration  scheme 

uJ+1  =  Lrf+(I-LrL)Uj. 

The  iteration  error  can  be  found  from  taking  G  -4  L1*  in  the  following  theorem,  which  is  well-known  in 
matrix  iterative  analysis. 

Theorem  3  Let 


uj+ 1  =  Gf  +  (/  -  GL)uj 

=  Uj+Gif-Luj)  (2.27) 

represent  an  iteration  scheme  for  uj  for  generic  G.  If  the  exact  error  in  a  .  is  given  by  ej  =  u  —  Uj  for 
u  =  L~l  f ,  then 

eJ+1  =  (I  -  GL)tj.  (2.28) 


«;+i  =  Gf+(I-GL)uj 

=  GL(L~1f)  +  uj  -  GLuj 
=  GLu  +  uj  —  GLuj 
—  GL(tt  —  v  j )  4-  u  j 
=  G  Lt  j  -T  u  j . 


Proof:  Equation  (2.27)  leads  to 
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so  that 


ej+ 1  =  «  “  «j+ 1 

=  u  —  i/?-  —  GLej 
—  Cj  —  GLe  j 
=  {I  -  GL)e-j. 

As  (/  —  L1  L)  is  idempotent,  its  norm  is  never  less  than  one,  and  the  iteration  scheme  does  not  improve 
the  error  with  each  iteration.  As  in  the  multigrid  case  (see  Lemma  2.1  of  [1,  p.  23]),  one  needs  to  combine 
L1  with  another  approximation  for  L~l .  With  the  conjugate  approximation  for  L'1  denoted  by  Z-1 ,  the 
iterative  scheme 

u0  -  0 

Mj+i  =  L~]  f  +  (I  -  L~l L)Uj 
uJ+1  =  L1*  f  +  (I  -  L1* L)uj+\ 

for  j  >  0,  for  which  the  error  satisfies 

ej+i  =  (I  -  L1* L)(I  -  L~1L)ej, 

converges  if  the  norm  of  (I  —  L~^  L)  is  sufficiently  small.  One  would  expect  that  if  Z^1  is  a  “good  enough” 
approximation  for  L~l ,  so  that  L~lL  is  “close  enough”  to  /,  then  this  norm  would  be  small.  If  Z-1,  is  a 
very  good  approximation  for  L~l ,  then  one  iteration  should  be  a  sufficient  approximation,  so  that  j  —  0  in 
(2.31)  leads  to 

L~l  «  L7*  +  (/  -  Z/*L)Z-1  (2.32) 

upon  taking  z/x  «  L”1  /  on  the  left-hand  side  and  dropping  /  (because  /  is  arbitrary).  Equation  (2.32),  or 
some  variant  of  it  based  on  (2.25)  or  (2.26),  can  form  the  basis  of  some  noniterative  model  reduction  hybrid 
methods. 

The  iteration  scheme  of  (2.30)  and  (2.31)  is  analogous  to  a  bilevel  multigrid  method  [1,  2],  with  (2.30) 
analogous  to  the  fine-grid-smoothing  step  and  (2.31)  analogous  to  the  coarse-grid-correction  step.  The 
operator  (I  —  L1  L )  is  analogous  to  the  coarse-grid-correction  matrix.  This  suggests  possible  variants 
of  (2.30)  and  (2.31),  such  as  using  multiple  iterates  of  (2.30)  both  before  and  after  the  (2.31)  iterate,  as 
is  usually  the  case  in  the  conventional  bigrid  method.  Even  at  this  abstract  level,  the  algebraic  essence 
of  the  bigrid  method  is  preserved  in  that  (2.12)  and  jV (I  -  l/*L)  =  TZ(Pd)  of  (2.15)  are  analogous  to 
A  [lfLhAh)  C  1Z(CG)  (implied)  and  N{CG)  =  7l{I}2h),  respectively,  of  Briggs  [2,  p.  79].  (Both  (2.13)  and 
(2.14)  are  satisfied  in  the  conventional  bi-grid  method,  as  will  be  seen  later.) 

2.3  Generic  Perturbational  Conjugate  Approximation 

A  generic  continuation  approach 

=  T(0-1|f  =  0_n  (2.33) 

for  T(  1 )  =  L  will  now  be  developed  as  the  conjugate  approximation  o f  L  ~ 1  for  use  in  (2.30)  or  (2.32).  As 
a  brief  synopsis,  in  continuation  methods  one  embeds  the  problem  to  be  solved,  denoted  generically  by 


(2.29) 

(2.30) 

(2.31) 


T(e)u  =  </(?). 
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into  a  continuum  of  problems 


T(()u  =  q(e). 


linked  by  a  parameter,  in  this  case  e.  The  problem  to  be  solved  is  placed  at  e  =  e,  where  t  is  usually 
normalized  so  that  7  =  1,  and  a  related  problem 


T(0)u  =  9(0). 

which  is  comparatively  easy  to  solve,  is  placed  at  (  —  0.  There  are  a  number  of  ways  in  which  the  computed 
continuation  from  e  =  0  to  e  =  I  for  the  generic  case  can  be  accomplished  [31, 53],  One  approach  is  to  carry 
out  a  perturbation  expansion  of  the  embedded  e  -problem  about  e  =  0,  carry  a  sufficient  number  of  terms  in 
the  expansion  for  accuracy,  and  evaluate  the  resulting  expansion  at  e  =  t.  This  is  an  old,  widely  used  version 
of  continuation  (Ref.  [54,  p.  245]),  and  it  is  the  version  that  is  used  in  this  report  (with  ?  =  1). 


If  the  t  expansion  of  T  ( e )  is  finite  for  generic  operator  T,  then  an  explicit,  generic  perturbational  expan¬ 
sion  for  T(£)_1  is  given  by  the  following  two  theorems. 

Theorem  4  Let  the  operator  T  have  the  finite  expansion 

./ 

T{()  =  E  (JTj  (2.34) 

j= o 

in  the  scalar  parameter  (  for  a  given  nonnegative  integer  J.  The  linear  component  operators  Tj  are  each 
assumed  to  be  independent  of  e,  and  To-1  is  assumed  to  exist.  Let  the  operators  and  i  L  each  have  the 
finite  expansions 

N 

r6(6)  =  JVrb,  (2.35) 

j= o 

for  a  given  nonnegative  integer  Ar,  b  — >  R  or  L,  where  the  component  operators  TRj  and  FL  j  are  each 
assumed  to  be  independent  of  e.  Define  the  right  component  operators  f by  the  recursion  relation 


TofTj  =  HNj  - 


E  HjkTknRj-k 


U=1 


and  the  left  component  operators  QLj  by  the  recursion  relation 


nLjT0  =  HXjrLj 

where  the  discrete  step  function  Hjk  is  defined  by 

H;k  =  { 


E  Hjk&j-kTk 

.7=1 


1  if  k  <  j 
0  if  k  >  j, 


(2.36) 


(2.37) 


(2.38) 


and  where  /  is  the  identity  operator  and  0  is  the  zero  operator.  The  operators 

n*(€)  =  !>'«*;• 

j=o 


(2.39) 
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and 


M 

nL{t)  =  Y.(JQL:i' 

j= o 


(2.40) 


for  a  given  nonnegative  integer  M,  satisfy 

Ma x{X.M+J) 

TQR  -rR  =  y 

j=M+ 1 


'E  HjkHM,_kTkQR^)  -  //v(I',?: 


L  \A— 0 


and 


respectively. 


nLr  -  rL  = 


Max(.\'j/+j) 


E 


E  HjkHM,  ,7,..  -  f/.v,!'1 

L \k= 0  / 


(2.41) 


(2.42) 


Appendix  B  provides  the  proof.  The  above  result  has  utility,  for  example,  when  the  right-hand  sides  of 
(2.41)  and  (2.42)  are  (4(f'1/+1).  In  this  case,  (2.41)  and  (2.42)  take  the  practical  forms 


T9r  -  Tr  =  0(eM+1 } 
QlT-Tl  =  0{eM+1 } 


(2.43) 

(2.44) 


It  is  assumed  that  even  if  the  perturbation  is  not  analytical,  the  results  still  have  value  in  terms  of  an  asymp¬ 
totic  series.  The  “of  the  order  of”  Landau  symbol  O  (Ref.  [8,  p.  8])  is  defined  by 


{  .4(6)  =  0(6*)  }  ►  {  limf_K, 


-4(c) 


fk 


<  OC 


} 


for  generic  A  of  a  normed  space  with  given  norm  1 1  1 1  and  a  given  nonnegative  integer  k.  Sufficient  conditions 
under  which  9L  and  9I!  coincide  are  provided  by  the  following  theorem.  (The  question  as  to  whether  these 
conditions  are  also  necessary  will  not  be  pursued  here.) 


Theorem  5  Under  the  hypothesis  of  Theorem  4  for  N  =  0,  assume  that  (2.36)  and  (2.37)  reduce  to 


nRm  =  H0mA  -  E  HmkATk9Rm-k 

k=  1 


and 


n1 


m 


=  H0mA  -  E  Hmk9Lm_kTkA. 


(2.45) 


(2.46) 


k=  1 


respectively,  for  some  operator  A  —  T0  1TR 0  =  TL0To  x.  For  each  m  >  0,  the  QRm  component  of  (2.45) 
and  the  corresponding  QLm  component  of  (2.46)  are  equal. 

The  proof  is  deferred  to  Appendix  C. 

2.3.1  Perturbational  Operator  Inversion 

An  important  special  case  of  the  previous  section  is  that  for  which  Ar  =  0  and  VR  =  TL  —  F0  =  /,  the 
identity  operator,  so  that  (2.36)  and  (2.37)  reduce  to 


nRm  =  H0mT0-1  -  E  HmkT(y  J/,0 

k=\ 


R 


(2.47) 
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nLm  =  H0mT0-1  -  Y,  HmkQLm-kTkT0-1 .  (2.48) 

k=i 

respectively,  and  (2.43)  and  (2.44)  together  imply 


Q  =  T(f)-1  +  C2(cA/+1).  (2.49) 

One  would  expect  this  to  be  the  case  for  T (e)  bounded-holomorphic  (Ref.  [55,  pp.  366, 419]),  for  example. 
The  superscripts  R  and  L  have  been  dropped  from  the  V  in  (2.49)  because  in  this  case  one  has  QL  =  QR, 
as  supported  by  the  following  corollary  to  Theorem  5. 

Corollary  2  Under  the  hypothesis  of  Theorem  5  for  TR  =  TL  =  I,  so  that  .4  =  Tq~  1 ,  assume  that  (2.41)  and 
(2.42)  imply  (2.43)  and  (2.44),  respectively.  For  each  j,  the  QRj  component  of  (2.47)  and  the  corresponding 
Ql  j  component  of  (2.48)  are  equal,  so  that  (2.49)  and 

M 

T(e)_1  =  Y  ejQJ  +  0{eM+l)  (2.50) 

j= o 

are  justified  for  Q  j  component  operators  recursively  defined  by  either  (2.47)  or  (2.48). 

Corollary  2  shows  that  the  components  given  by  (2.47)  can  be  thought  of  as  those  of  the  e-expansion  of 
J"1 .  The  only  inversion  used  in  computing  T(e)”1  for  all  6  is  that  for  To. 

The  results  (2.47)  can  be  tested  against  known  results  for  generic  perturbational  operator  inversions. 
The  .7  =  1  subcase  of  (2.50),  for  which  (2.47)  reduces  to 

Qo  =  T0 -1  (2.51) 

=  [-T0-1Ti]Om_i  for  m  >  1,  (2.52) 

is  equivalent  to  the  well-known  Neumann  series  (Ref.  [55,  pp.  30, 32]).  To  see  this,  take  J  =  1  in  (2.34)  to 
get 


T  =  To  +  d\ 

=  To  [/  +  eTcT1!)] 

=  To  [/-€(— To-'Tj)], 


so  that 

T"1  =  [/-e(-To"1ri)]"1T0-1.  (2.53) 

The  Neumann  series  results  by  the  use  of  the  “binomial”  operator  expansion 

CO 

[/-  e.4]-1  =  Y  (2-54> 

?n= 0 


for  bounded  linear  operators  .4  (Ref.  [56,  p.  375]).  (The  expansion  (2.54)  is  convergent  when  the  magni¬ 
tude  of  c  is  less  than  the  inverse  of  the  norm  of  A ,  assuming  a  Banach  space  setting.)  Using  (2.54)  with 
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A  =  —To  lT\  in  (2.53)  leads  to 


T-i  = 


| -To 

l  m  =0 

E  <”7  [— To_1Ti 

'77=0 

E 


‘T, 


To'1 


To” 


m=0 


where 


fim  -  [-ro-'T! 

This  expression  for  f2m  is  equivalent  to  (2.51)  and  (2.52). 
For  the  J  =  2  subcase  of  (2.50),  (2.47)  leads  to 


To 


-l 


fto  =  To"1 
fii  =  -T0"  i  i-t0 
Qm  =  —To  1[Txf2m_i  +  T2fim_2] 


'“TrTr1 


(2.55) 


(2.56) 


(2.57) 

(2.58) 

(2.59) 


for  zn  >  2.  This  agrees  with  the  first  few  terms  for  T  1  given  by  Kato  in  (6.16)  of  (Ref.  [55,  pp.  420]). 
Note  that  (2.57)  through  (2.59)  reduce  to  (2.51)  and  (2.52)  for  T2  =  0,  as  it  should. 


232  Outline  of  Procedure  for  Applying  Formalism 

To  complete  the  development  of  the  generic  continuation  (2.33)  for  T ( 1 )  —  L  as  the  conjugate  approxima¬ 
tion  of  L~l ,  take  e  — y  1  in  (2.50)  to  get 


M 

r(€)-1|f  =  0_>1  =  E«;-  (2-6°) 

3= 0 

As  (2.47)  or  (2.48)  are  equally  valid  by  Corollary  2,  (2.47)  will  be  used  for  definiteness.  With  the  superscript 
R  dropped  from  QR  in  (2.47)  because  QL  =  the  Qf  s  in  (2.60)  are  recursively  given  by 

J 

Qm  —  H0mT0  1  -  ^2  Ttmk'To  1TkQm-k-  (2.61) 

h= i 

The  steps  to  implement  the  formalism  of  this  paper  using  this  generic  perturbation  as  the  conjugate  approx¬ 
imation  are  summarized  as: 

•  Construct  the  linear  idempotent  operators  Pr  and  Pcj  such  that  the  relevant  stimuli  and  response  sub¬ 
sets  are  given  by  1 Z (Pr)  C  Tv ( L )  and  7 Z(Pd)  C  V  ( L ) ,  respectively,  where  L  is  the  system  operator  for 
the  system  model  (2.1).  It  is  advantageous,  but  more  work,  to  fix  one  of  either  Pr  or  Pfj  and  determine 
the  other  from  either  (2.21)  or  (2.23). 

•  For  Lfjj  given  by  (2.10),  determine  Lcj j1  such  that  it  satisfies  either  (2.13),  (2.14),  or  both.  Construct 
Lr  from  this  LC:jfl  using  (2.11). 


•  Embed  L  in  a  a -continuation  T(e)  such  that 
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-  To-1  exists 

-  To-1  is  much  more  easily  computed  than  L~l 

-  T(l)  =  L,  and 

-  T  has  the  expansion  (2.34)  for  some  J. 

•  The  continuation  should  complement  the  model  reduction  approximation  L1  in  that  it  should  lead  to 
a  reasonably  accurate  approximation  to  L~l  outside  of  lZ(Pr)  and  7 Z(Pd). 

•  The  conjugate  approximation  L"1  is  determined  by  (2.33)  and  (2.60)  for  Q:I ’s  recursively  given  by 
(2.61)  for  some  chosen  value  of  M. 

•  The  solution  u  to  (2.1)  is  constructed  either  from  u  =  L~l  f  with  L~l  approximated  by  (2.32),  or 
iteratively  using  (2.29)  through  (2.31)  (or  some  variant).  There  will  probably  be  a  tradeoff  between 
M  and  the  number  of  inner  iterations  of  (2.30);  accuracy  may  require  many  iterations  of  (2.30)  per 
evaluation  of  (2.31)  for  small  A/,  or  only  one  iteration  of  (2.30)  per  evaluation  of  (2.31)  for  large 
enough  M. 


The  next  section  shows  that  there  are  several  important  classes  of  model  reduction  methods  for  which  the 
second  step,  in  which  L?  f  j1  is  determined  to  satisfy  either  (2.13),  (2.14),  or  both,  is  already  “worked  out.” 
The  operator  L€j j1  is  constructed  from  where  Lre({  is  a  “reduced”  version  of  L.  The  construction 

of  Le  f  f 1  from  Lrec/_1  and  the  specific  definition  of  Lred  varies  from  one  class  of  model  reduction  method 
to  another.  For  these  cases,  the  bulk  of  the  effort  in  the  second  step  above  reduces  to  the  computational  work 
of  obtaining/applying  £red_1 . 

3.  CLASSES  OF  METHODS  ENCOMPASSED 

The  formalism  of  this  report  encompasses  the  reduced  basis,  substructuring,  and  smoothing/homogenization 
methods  of  model  reduction.  To  maintain  as  general  a  setting  as  possible  for  the  development  of  the  reduced 
basis  and  substructuring  methods,  L  of  (2.1)  is  taken  to  be  a  square  system  matrix  with  operator  compo¬ 
nents,  as  in  Ref.  [20].  In  keeping  with  this,  the  definition  of  the  multiplication  of  two  arbitrary  matrices  is 
generalized  to 

{AB)a  =  (2U) 

k 

for  compatible  matrices  .4  and  B ,  where  the  symbol  o  denotes  mapping  composition,  that  is,  AikoBkj 
applied  to  some  function  g  is  interpreted  as  Aik{Bj~j{g)).  The  associative  and  distributive  laws  of  ordinary 
matrix  algebra  carry  over  to  this  more  general  setting  for  the  case  of  linear  operator  components.  For  L  to 
preserve  compatibility  with  multiplication  by  ordinary  matrices,  and  to  associate  each  component  of  u  with 
one  degree-of-freedom,  the  domain  and  range  of  each  of  the  operator  components  consist  of  scalar- valued 
functions.  In  the  special  case  where  B  in  (3.1)  is  a  matrix  of  such  functions,  the  o  in  (3.1)  is  interpreted 
to  mean  Aj^oB^j  =  Ai^(B^}).  Ordinary  matrices  are  special  cases  for  which  each  component  operator 
consists  of  scalar  multiplication  by  a  fixed  scalar  value.  Viewed  as  a  single  matrix,  the  L  of 

L-K  +  4  +  ilW  <3'2) 

is  a  less  trivial  example  of  an  operator-component  matrix,  each  component  being  an  ordinary  differential 
equation  operator.  In  mechanics  the  A/,  C1 ,  and  K  denote  the  mass,  damping,  and  stiffness  matrices,  respec¬ 
tively. 
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3.1  Reduced  Basis  Methods 

Reduced  basis  methods  consist  of  an  approximation  to  the  original  governing  equations  over  subspaces, 
each  of  which  is  spanned  by  a  small  number  of  basis  vectors.  Let  the  columns  of  the  matrix  $d  consist  of 
a  set  of  basis  vectors  for  a  chosen  subspace  of  the  domain  of  L.  The  span  of  the  columns  of  represents 
the  subspace  of  relevant  system  responses  u  for  (2.1).  Similarly,  let  the  columns  of  the  matrix  <3>r  consist 
of  a  set  of  basis  vectors  for  a  chosen  subspace  of  the  range  of  L.  The  span  of  the  columns  of  3>r  represents 
the  subspace  of  relevant  system  stimuli  /  for  (2.1).  The  matrices  <FC/  and  <E»r  are  of  the  same  dimensions. 
The  accuracy  of  any  particular  reduced  basis  method  is  largely  dependent  on  the  particular  choice  of  basis 
in  and  such  choice  largely  distinguishing  between  specific  methods.  Let  Tr  and  Td  be  two  square 
matrices  whose  size  is  equal  to  the  number  of  rows  of  and  <I>r,  such  that  $r^rr$r  and  $d^d$d  are  both 
nonsingular.  The  superscript  |  denotes  the  adjoint  (transpose  with  complex  conjugation).  Take  Pr  and  Pd 
as 


for 


pr  =  rvr/  0.3) 

Pd  =  r/rd  (3.4) 

r,7  =  $r[$rtrr$r]-1$rt  (3.5) 

r  /  =  ^[^rA]-1^/.  (3.6) 


Note  that  IT,.7  and  I^7  as  defined  by  (3.5)  and  (3.6)  are  outer  generalized  inverses  of  I\  and  Tj,  respectively. 
By  Theorem  1 ,  Pr  and  Pd  are  idempotent  as  required  by  Theorem  2.  It  is  also  readily  shown  that 


$rf  Pr  =  $rf  (3.7) 

Pd*d  =  (3.8) 


The  Lej j  defined  by  (2.10)  becomes 

Uss  =  rr$,[$,.trr$,.]-1Lrerf[$/r^rf]-1$/rd 

in  this  case,  where  the  reduced  version  of  L  for  this  class  of  methods  is  defined  to  be 

Lred  —  &r^L$d- 


(3.9) 


(3.10) 


The  size  of  Lred  is  equal  to  the  number  of  reduced  basis  vectors,  that  is,  the  number  of  columns  of  <hr  (or 
$rf).  The  essence  of  the  reduced  basis  method  is  to  use  Lred~l  instead  L~l  to  solve  (2.1).  Therefore,  the 
number  of  such  basis  vectors  must  be  kept  small  compared  to  the  original  number  of  DOF  for  the  method 
to  be  efficient.  The  number  of  <3>r ’s  (or  <L/’s)  columns  must  be  significantly  less  than  the  number  of  its  rows 
if  worthwhile  problem  size  reduction  is  to  occur. 


It  is  readily  verified  that  Le  /  f1  given  by 

Lf//  =  $,,Lrerf-1$rt  (3.11) 

is  an  outer  generalized  inverse  of  the  Lejj  of  (3.9)  that  also  satisfies 

Ltff  Lfjj  = 

LtJfL'ff1  = 


Pd 

Pr- 


(3.12) 

(3.13) 
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Since  Pr  and  Pi  are  idempotent,  (3.12)  and  (3.13)  imply  (2.13)  and  (2.14),  respectively.  It  is  also  readily 
verified  that  iJ”  =  Lej /  for  Lr  given  by  (2.1 1),  so  that 

Lr  =  (3.14) 

when  (3.10)  and  (3.11)  are  used.  The  classic  Rayleigh-Ritz  reduced  basis  approximation  consists  of  (3.14) 
for  the  special  case  of  4>(/  =  4>,.  =  4>,  appropriate  for  self-adjoint  L. 

3.1.1  Modal  Method  for  Linear  Mechanics  FEM 

As  a  common  and  important  application  of  reduced  basis  methods,  consider  the  linear  dynamic  mechanics 
class  of  finite-element  models  consisting  of  second-order  ordinary  differential  equations  with  real,  constant 
symmetric  matrix  coefficients.  The  L  of  (2.1)  becomes  (3.2),  for  which  the  typical  M  and  K  are  at  least 
positive,  if  not  positive  definite.  The  modal  reduced  basis  method  as  typically  applied  to  such  systems  can  be 
derived  as  the  subcase  for  which  T,  =  Td  =  M  and  4>d  =  4>r  =  4>,  where  the  columns  of  4>  are  normalized 
such  that 

4>TAf4>  =  /.  (3.15) 

The  superscript  T  denotes  the  transpose,  the  matrices  being  real.  For  this  case,  the  matrices  Pr  and  Pd 
reduce  to 

Pr  =  M  4>4>T 
Pd  =  4>4>rA7 

The  modal  method  is  usually  accompanied  by 

•  the  constraint  that  4>TL4>  be  diagonal,  and 

•  the  constraint  (2.23),  so  that  (2.24)  is  true.  This  means  that  /  =  Prf  =  L(Lr  f)  for  all  /  in  the  range 
of  Pr,  for  which  u  =  1. '  /  is  an  exact  solution  to  (2.1). 

Substituting  (3.16)  and  (3.17)  into  (2.23),  matrix  multiplying  the  result  on  the  right  by  4>,  and  then  using 
(3.15)  on  the  result  gives 

(/  -  A/$4>T)A4>  =  0. 

This  leads  to 

A’4>  =  A/4>T2 

C'4>  =  A/4>J, 

where 

T2  =  4>rA'4>  (3.20) 

3  =  (3.21) 

for  K  and  C,  (I  -  A/4>4>T)A/4>  =  0  being  an  identity  for  A/  by  (3.15).  The  constraint  that  4>TL4>  be 
diagonal  reduces  to  the  constraint  that  Y2  and  /3  be  diagonal. 

If  $  is  determined  so  as  to  satisfy  (3.18),  then  the  columns  of  $  are  eigenvectors  satisfying  the  real 
eigenproblem  (3.18)  for  which  the  eigenvalues  lie  along  the  diagonal  of  the  diagonal  matrix  Y2.  The  eigen¬ 
vectors  are  mass  normalized  by  (3.15).  (One  usually  takes  4>  such  that  I\  and  A/  are  positive  definite  over 


(3.18) 

(3.19) 


(3.16) 

(3.17) 
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the  span  of  the  columns  of  <J>,  so  that  the  I  in  (3.15)  and  the  square  of  T2  in  (3.20)  both  make  sense.)  The 
relation  (3.19)  is  a  constraint  on  C.  It  can  be  put  into  a  more  conventional  form  by  the  following  argument. 
Substituting  (3.17)  in  CPd  leads  to  CPd  —  M<&  M>7 .1/  upon  using  (3.19).  Substituting  this  into  the  identity 
C  —  C'Pd  +  C (/  —  Pd)  produces  the  constraint 

C  =  (il/$)/J($rA/)  +  C(I  -  Pd)  (3.22) 

on  C .  Take  (3.21)  in  (3.19),  transpose  the  result,  and  then  use  CT  =  C\  MT  =  M9  and  (3.17)  to  get 
<&TC  —  $TC($$T M)  =  $TCPd,  so  that  <1>TC(/  -  Pd)  =  0.  Matrix  multiplying  this  on  the  left  by  3/$ 
gives  PrC{I  -  Pd)  =  0  by  (3.16),  so  that  the  identity  CJ(I  -  Pd)  —  PrC(I  -  Pd)  +  (I  -  Pr)C(I  -  Pd) 
gives  C{I  -  Pd)  —  (/  -  P,  )C(I  -  Pd)  =  [I  -  M$$r)C(/  -  $<hTM).  Substituting  this  into  (3.22)  and 
using  3/t  =  il/  finally  leads  to  the  general  constraint 

C  =  {AI$)i3{M<S>)t  +  (/  -  M$$t)C(I  -  M$<S>t)t  (3.23) 

on  the  damping  matrix.  The  constraint  (3.19)  implies  that  (3.23)  and  the  converse  is  true  as  well.  If  M 
is  nonsingular  and  one  has  a  complete  set  of  modes,  $  then  being  square,  then  one  would  expect  (I  - 
MQ$t)  =  0  in  (3.23)  andC  =  (M$)t3{M$)T.  In  practice,  the  (/  -  M$$T)C(I  -  M$$T)T  component 
of  C  is  usually  neglected,  even  if  one  has  an  incomplete  set  of  modes,  so  that  C  ~  (M<fr)i3(M&)T .  Also, 
the  3  matrix  is  usually  specified.  A  common  choice  for  /3  in  such  cases  is 

8  =  2CT,  (3.24) 

where  (  is  a  diagonal  matrix  with  damping  ratio  values  along  its  diagonal.  Under  the  condition  (3.15)  with  a 
complete  set  of  modes,  Wilson  and  Penzien  (see  (17)  through  (19)  of  Ref.  [57])  obtain  C  = 
with  3  given  by  (3.24).  This  constraint  on  C  generalizes  Rayleigh  and  proportional  damping. 


3.1.2  The  Conventional  Biscale  Case  of  Multigrid 

In  the  biscale  case  of  multigrid  methods,  the  coarse-grid  correction  is  essentially  a  reduced  basis  method. 
Using  the  notation 


L  — ^  f  /  jj 
u  — y  u  /, 

i  -»>  fh 

of  Ref.  [1,  pp.  18-28],  (2.1)  represents  the  fine  grid  problem.  The  bigrid  coarse-grid  correction  of  Ref.  [1, 
p.  21]  is 


u 


j+i 

h 


<  +  LlX[f  ~  Lhu{ ] 
L1* f  +  (/  -  L1* Lh)a{ 


for  L1*  given  by 


I^Lh-1!? 


h  • 


This  coarse-grid  correction  is  of  the  form  (2.31).  The  above  expression  for  Lr  implies 


LI"  =  Il}IVl?LhIhH)-1lf? 
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when  using 


Lh  =  J\ 


H  r  rh 
Lh*H 


from  (2.21)  of  Ref.  [1,  p.  27].  This  is  in  the  form  of  (3.14)  upon  taking 


*d  =  Ih 

=  it 


The  iff  and  Pjq  are  the  restriction  and  prolongation  operators,  respectively,  required  for  intergrid  informa¬ 
tion  transfer.  The  matrix-splitting  class  of  error  smoothers  used  in  multigrid  can  be  recast  as  continuation 
methods.  As  an  example,  splitting  Lh  =  A  -  U,  where  A  is  lower  triangular  and  U  is  upper  triangular  with 
zero  diagonal,  suggests  the  continuation  T(e)  =  A  —  eU  with  T(l)  =  Lh-  For  this  case,  the  zeroth-order 
(e  — /  0)  version  of  Z-1  in  (2.33)  leads  to  the  well-known  Gauss-Seidel  matrix  iterative  scheme  for  (2.30). 


3.2  Substructuring  Methods 

In  substructuring  methods,  one  tries  to  approximate  the  response  of  (2.1)  for  a  chosen  subset  of  the  DOF 
of  (2.1);  this  subset  being  denoted  as  master  DOF.  The  remaining  DOF  are  correspondingly  denoted  as 
the  slave  DOF.  For  many  substructuring  applications,  the  slave  DOF  are  associated  with  a  collection  of 
substructures  that  are  coupled  only  through  the  master  DOF,  the  master  DOF  being  associated  with  a  main 
coupling  structure.  Permuting  the  DOF  into  master  and  slave  subsets  can  be  represented  mathematically  by 
a  (real,  orthogonal)  permutation  matrix  P,  where 

pT  =  p~l  (3.25) 


and  the  superscript  T  denotes  the  transpose.  The  permutation  matrix  V  is  defined  to  permute  u  of  (2.1),  so 
that  the  resulting  first  M  components  of  V  u ,  collectively  denoted  by  um ,  are  associated  with  the  M  master 
DOF,  and  the  remaining  components,  collectively  denoted  by  us,  are  associated  with  the  slave  DOF.  The 
convention  adopted  in  this  report  is  that  V,  by  definition,  gathers  the  master  DOF  into  the  upper  part  of  V  u , 
so  that 


Pu  = 


U-m 

Us 


(3.26) 


and  similarly  for  any  other  column  matrix  of  the  same  size  as  u.  Using  this  convention,  a  choice  of  sub¬ 
structuring  is  completely  determined  mathematically  by  the  permutation  matrix  V  and  the  number  of  master 
DOF  M .  This  can  be  extended  to  square  matrices  of  the  size  of  u  as  well,  and  in  particular 


L  =  VLV~l. 

for  L  leads  to  the  block  partitioned  form 


(3.27) 


(3.28) 


That  the  L'P~l  part  of  (3.27)  permutes  the  columns  of  L  can  be  seen  by  LP~ 1  =  LP1  =  [PLT]T ,  where 
P Lt  permutes  the  rows  of  LT .  The  Lmm  of  (3.28)  is  a  M  x  M  matrix  consisting  only  of  those  components 
of  L  relating  master  DOF  to  master  DOF.  A  similar  statement  is  true  for  Lss  with  respect  to  the  slave  DOF. 
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From  the  point  of  view  of  this  report,  substructuring  methods  are  also  reduced  basis  methods  but  with 
a  different  structure  from  that  of  Section  3.1.  In  both  cases,  the  span  of  the  columns  of  <J>  determines  the 
subset  of  'JZ(L)  of  ‘relevance.”  In  this  case,  however,  the  projector  Pr  satisfies 


PA>r  =  <lv 


(3.29) 


for  some  choice  of  as  opposed  to  (3.7).  In  addition,  the  admissible  <I\.  ’s  from  which  one  may  choose  are 
restricted  to  a  class  which  is,  in  some  sense,  compatible  with  one’s  choice  of  substructuring  for  the  system. 
This  notion  is  made  more  precise  by  the  following  definition. 

Definition  2  For  a  given  number  of  master  DOF  M  and  a  given  choice  of  the  permutation  matrix  V,  let 
the  columns  of  <3>r  span  a  subset  of  7 Z(L).  The  reduced  basis  consisting  of  the  columns  of  <I>,  is  said  to  be 
compatible  with  the  substructuring  M  and  V  if 


t 


det ([$,.], n  [<frr]TO)  #  0 


for  the  M  x  K  matrix  1  <  K  <  M ,  where 

*r  = 

for  given  by 


r]n 

[*r]s 


=  T$r. 

The  subscripts  m  and  s  denote  master  and  slave,  respectively. 


(3.30) 


(3.31) 


(3.32) 


The  following  theorem  gives  this  reduced-basis-compatible  version  of  substructuring  in  the  context  of  the 
formalism  of  this  report. 

Theorem  6  Let  the  reduced  basis  consisting  of  the  columns  of  <I>,  be  compatible  with  the  substructuring 
M  and  V  as  in  Definition  2.  Let 


n-y  = 


0 


(3.33) 


represent  shorthand  notation  for  a  matrix  function  of  its  submatrix  7  such  that  square,  idempotent  II-  is  the 
same  size  as  L.  The  Imm  of  (3.33)  is  an  ,V(  x  M  identity  matrix.  If  the  matrix  a  is  defined  by 


Q  =  [$r]a  ([#,-]„,  t[$,]m)“1 

from  (3.31)  and  (3.32),  then  the  idempotent  matrix  Pr,  defined  by 

pr  =  v~ln0v. 

satisfies  (3.29).  If  3  satisfies 

[■^ss  ^  —  [oZ/m m  Lsm], 

then  (2.13),  (2.14),  and  (2.23)  are  satisfied  for  Pj,  Lre,i,  and  Lf  j j 1  given  by 

Pd  =  V~%)V. 

bred  ”  Lmm  T  Lms3 


L, 


JJ 


=  V~lQLP. 


(3.34) 

(3.35) 

(3.36) 


(3.37) 

(3.38) 

(3.39) 
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where 


so  that 


0  = 


0'  = 


Lred  0 

0  0 


red 


-l 


In  this  case,  Lrfj  is  the  reduced  version  of  L  and  Lr  of  (2.11)  and  becomes 

L1*  =  p-'iuVP. 


(3.40) 


(3.41) 


(3.42) 


Appendix  D  contains  the  proof.  The  hypothesis  of  Theorem  6,  if  satisfied,  implies  (2.14)  and  (2.23),  which 
in  turn,  by  Theorem  2,  implies  (2.24).  As  shown  previously  (just  after  Theorem  2),  (2.24)  implies  that 
L1*  f  is  an  exact  solution  to  (2.1)  for  all  /  such  that  Prf  =  f.  The  special  case  of  a  =  0  in  (3.35)  corre¬ 
sponds  to  conventional  substructuring  and  its  variants  as  currently  practiced,  with  the  zero  frequency  limit 
corresponding  to  Guyan  condensation.  For  o  =  0,  Prf  =  /  implies  /  =  P-1lloP/,  so  that 


f  =  V~ 1 


fm 

0 


(3.43) 


The  slave  DOF  hence  cannot  be  loaded  for  conventional  substructuring  if  L1"  f  is  to  retain  its  exact  solution 
status  for  that  case.  The  idea  of  using  a  as  in  (3.35)  and  (3.36)  has  been  presented  before  [52,  20,  49],  but 
its  use  in  connecting  reduced  basis  methods  with  substructuring  methods  is  emphasized  here.  There  is  some 
similarity  in  the  use  of  (3.34)  for  a  and  the  Modal  Reduction  method  (see  (10)  of  Ref.  [21,  p.  327],  for 
example),  which  only  considers  the  case  in  which  the  reduced  basis  consists  of  system  modes.  The  Modal 
Reduction  method,  however,  builds  the  reduced  system  4>~L<P  from  the  transformation  matrix  (7  aT)T,  and 
the  similarity  ends  with  (3.34). 


3.2.1  Reduced-Basis/Substructuring  Relationship 

If  the  number  of  columns  of  <f\,  is  equal  to  M ,  the  number  of  master  DOF,  then  the  reduced  basis  approach 
of  Section  3.1  and  the  version  of  substructuring  given  by  Theorem  6  can  be  directly  related  to  each  other. 

Theorem  7  Assume  the  hypothesis  of  Theorem  6.  Take  the  columns  of 

=  V~lTla^  (3.44) 

as  the  reduced  basis  for  1Z{L),  and  similarly,  take  the  columns  of 

=  p-1n/3(  j  (3.45) 

as  the  reduced  basis  for  V(L),  where  3  satisfies  (3.36),  and  [$,.]„,  and  [$</]m  are  each  square  and  nonsin¬ 
gular.  For  comparison  purposes,  denote  the  substructuring  version  of  L 1  from  (3.42)  as  Z;*u6.  For  <I\  and 
given  by  (3.44)  and  (3.45),  respectively,  denote  the  reduced  basis  version  of  L1'  from  (3.14)  as  Lr'rh. 
The  i 1  and  L,xsuh  are  directly  related  to  each  other  by 


(3.46) 
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The  proof  is  deferred  to  Appendix  D.  Note  that  of  (3.46)  is  idempotent.  This  im¬ 

plies  that  L1*!,  and  coincide  over  the  range  of  <f>r($rt4>r)_14>rt.  Applying  (3.46)  to  <f>r  and  using 

to  the  result  leads  to  This  implies  that  L7*6  and  Z,7*ub 

coincide  over  the  span  of  the  columns  of  <Fr. 

3.3  Smoothing/Homogenization  Methods 

Fishman  and  McCoy  (Ref.  [13,  pp.  47-48])  have  unified  smoothing  and  homogenization  under  one  formal¬ 
ism.  Proof  that  the  Fishman  and  McCoy  formalism  is  a  subformalism  of  this  report  constitutes  proof  that 
smoothing  and  homogenization  are  encompassed  by  this  report.  Smoothing  methods  presuppose  the  exis¬ 
tence  of  a  “comparison  operator”  Lq,  a  user-chosen  approximation  to  the  system  operator  L  of  (2.1).  The 
term  is  taken  from  applications  involving  the  smoothing  of  heterogeneous  material  response  in  which  the 
L  and  Lq  typically  represent  the  system  operator  of  a  heterogeneous  media  and  an  associated  “comparison 
media,”  respectively.  The  comparison  media  is  often  assumed  to  be  spatially  homogeneous  with  constant 
constitutive  parameters  which  are  “close  enough”  to  the  spatially  fluctuating  ones  of  L  to  make  the  differ¬ 
ence  between  L  and  Lq  a  perturbation.  In  the  abstract  case,  the  comparison  operator  Lq  is  a  linear  operator 
that  is  presumed  to  satisfy,  by  definition,  the  conditions 

M{Lo)  =  {0}  (3.47) 

L0P  =  PrL0  (3.48) 

for  Pr  of  Theorem  2  and  P,  a  linear  idempotent  operator  for  which  7 Z(P)  C  V(L).  Condition  (3.47) 
guarantees  that  L0~l  exists.  In  applications  involving  the  smoothing  of  the  linear  response  of  stochas¬ 
tic  heterogeneous  materials,  for  example,  the  P  and  Pr  are  usually  both  taken  to  be  equal  to  a  common 
ensemble-averaging  projector.  A  general  formalism  for  smoothing  is  summarized  in  the  following  theorem. 


Theorem  8  Let  Lq  be  a  linear  operator  satisfying  the  conditions  (3.47)  and  (3.48).  Define  SL  =  L  -  L0  for 
L  of  (2.1),  so  that 


If  L  1  exists  for 


L  —  Lq  T  SL. 


(3.49) 


and  Pc{  is  defined  by 
then  (2.23)  and 


Z=L0  +  {I-Pr)6L 

(3.50) 

Pd  =  [I  -Z-'il  -  Pr)6L]P. 

(3.51) 

PPd  =  P 

(3.52) 

to 

II 

i3 

(3.53) 

Leff  =  {l0  +  Pr(SL)[T  -S-1!/  —  Pr)(SL)]}p 

=  P,.{l0  +  [I  -  (6L)^-l(I  -  Pr)](6L)p}. 

(3.54) 

for  Lefj  given  by  (2.10). 


Corollary  3  Under  the  hypothesis  of  Theorem  8,  if  the  constraint 


Pr(SL)P  =  0 


(3.55) 
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is  satisfied,  then 


Pd 


L<// 


[I  -  E-18L]P 
Lo  -  Pr{8L)E~1(8L) 


P 


Pr  Lo  -  (SL)E~1(8L)P 


(3.56) 

(3.57) 


for  Lf  jj  given  by  (2.10). 

Proof:  Let  G  -  1  ( I  - Pr )  so  that  EG  =  (/- Pr ) ,  and  0  =  Pr  ( /- Pr )  =  Pr EG  =  Pr  [70+(7  -  Pr ) 6L]G  = 

L0PG  by  (3.50)  and  (3.48).  This  means  that  7 Z(PG)  C  A r(L0),  so  that  PG  =  0  by  (3.47).  This  proves 

PE-1  {I-  Pr)  =  0.  (3.58) 


Using  this  and  (3.51)  in  PPd  =  P[I  -  E~1(I  -  Pr)8L]P  =  P  proves  (3.52).  Using  (3.51)  and  (3.52)  in 
Pd 2  =  [I  -  E~1{I  -  Pr)SL]PPd  =  [I  -  S_1(/  -  Pr)6L]P  =  Pd  proves  (3.53).  To  prove  (2.23),  take 
u  €  TZ(Pd)  so  that  a  =  Pju  =  [I  -  E~1(I  -  Pr)SL]Pu  by  (3.51)  and  (3.53).  (An  idempotent  operator 
acts  as  the  identity  over  its  range.)  This  leads  to  (/  —  P)u  =  Pr)(8L)Pu,  so  that  E(I  — 

P)u  =  -(/  -  Pr)(5L)Pu.  However,  (/  -  Pr)L(I  —  P)  =  (I  -  Pr)[L0  +  SL]{I  -  P)  =  L0{I  -  P)  + 
(. I  -  P,-)8L(I  -P)  =  [Lo  +  (/  -  Pr)SL](I  -P)  =  E(I  -  P)  by  (3.49),  (I  -  Pr)L0  =  L0(I  -  P)  from 
(3.48),  the  idempotent  property  of  P  (and  of  I  -  P),  and  (3.50).  Substituting  this  into  the  previous  result 
gives  (I  -  P,  )L(I  -  P)u  =  -(/  -  Pr)(8L)Pu.  The  interim  results  P,-L(I  -  P)  =  Pr{8L)(I  -  P)  and 
(/  -  Pr)LP  =  (I  -  Pr)(6L)P  follow  from  P,  Lq(I  -  P)  =  LoP(I  —  P)  =  0  and  (/  —  Pr)LoP  —  {I  - 
Pr)PrLo  =  0,  respectively.  Substituting  the  latter  interim  result  into  the  previous  result  gives  (I  —  Pr)L(I  - 
P)u  =  -(/  -  Pr)LPu,  so  that  0  =  (/  -  Pr)L(I  -  P)u  +  (/  -  Pr)LPu  =  [I  -  Pr)Lu.  Since  u  €  K{Pd) 
was  arbitrary,  (2.23)  is  proven.  Finally,  using  (3.51)  in  (2.10)  for  Lef  j  leads  to 

Lf//  =  PrL[I  -E-X(I  -  Pr)8L\P 

=  PrLP  -  PrLE~1(I  —  Pr)(8L)P 
=  Pr[L0  +  SL]P  -  PrL(I  -  P)S-1  (7  -  Pr)(SL)P 
=  LoP  +  Pr(8L)P  -  Pr(8L){I  -  P)E-\I  -  Pr)[8L)P 
=  L0P  +  Pr(SL)P-Pr(8L)E-1(I-Pr)(8L)P  (3.59) 

when  using  (3.49),  (7 -  P)^1  (7  -  Pr )  =  S"1  (7  -  Pr)  from  (3.58),  (3.48),  PrL(I -  P)  =  Pr  (47) (7 -  P), 
and  (7  -  P)L_1  (7  —  Pr)  =  E~1(I  -  Pr)  again,  respectively.  The  two  right-hand  sides  of  (3.54)  are  just 
rearrangements  of  (3.59),  where  (3.48)  was  used  in  the  latter. 


The  Lfj j  of  Theorem  8  takes  either  of  the  forms 

Lcfj  =  [Lred)dP 

=  P(Lrcd)r  (3-60) 

when  defining  the  reduced  system  operators 

(Lrtd)d  =  Lo  +  P(SL)[I-E-1(I-Pr)(6L)]  (3.61) 

{Lrfd)r  =  7o  +  [7  —  (8L)E~1{I  —  Pr)]{8L)P .  (3.62) 

=  L0-  Pr(8L)E~1(8L)  (3.63) 

=  L0-  {8L)E~1(8L)P  (3.64) 


which  simplify  to 


(Lr"i)d 

red)  r 
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when  (3.55)  is  satisfied.  Define  (Z,e//7)r  and  as 

(Lcjj1),.  =  (Lred);lPr  (3.65) 

(Pf/)d  =  P{Lr<id)dl-  (3.66) 

so  that  the  choice  of  ( Lre,i)r  for  the  reduced  operator  leads  to 

Lf.jj{L(jjI)r  =  Pr,  (3.67) 

which  gives  (2.14)  when  operating  on  (3.67)  from  the  right  with  Pr.  The  choice  (Lre(})d  leads  to 

(Pff^dPfj  =  P  (3.68) 

which  gives  (2.13)  when  operating  on  (3.68)  from  the  left  with  Pd.  Operating  on  (3.67)  on  the  left  by  (3.65) 
shows  that  [Ltj /),.  is  an  outer  generalized  inverse  of  Lfj /.  Similarly,  operating  on  (3.68)  on  the  right  by 
(3.66)  shows  that  (Lej  fT)d  is  an  outer  generalized  inverse  of  Ltj  j.  Taking  L1*  from  (2.11)  leads  to 

L r  =  PdLred~lPr  (3.69) 

for  generic  Lred;  that  is,  either  ( Lred)r  or  ( Lred)d  can  be  substituted  for  Lred  in  (3.69). 


If  /  satisfies  the  constraint 
so  that  /  €  1Z  ( Pr ) ,  then 


Prf  =  /• 

(3.70) 

u  =  Lrf 

(3.71) 

is  possibly  an  exact  solution  to  (2.1).  For  the  choice  of  (Lred)r  as  the  reduced  operator,  (2.23)  from  Theo¬ 
rem  8  and  (2.14)  together  imply  (2.24)  by  Theorem  2.  Assuming  that  {Lrtd)~l  f  exists,  (2.24)  and  (3.70) 
imply  that  (3 .7 1 )  is  a  solution  to  (2.1).  For  the  choice  of  ( Lred )  d  as  the  reduced  operator,  let  g  be  the  solution 
to 

{Lred)d9  =  P.f  (3.72) 


subject  to  the  constraint  that 


P(J  =  fir 


(3.73) 


so  that,  if  it  exists,  g  €  P{{Lred)d  1Pr)n72(P).  The  problem  (3.72)  is  the  “reduced  version”  of  (2.1)  for  this 
case.  Define  u  by 


u  =  Pig- 


(3.74) 


so  that  (3.71)  follows  from  (3.69)  with  (Lred)d  substituted  for  Lred.  Operating  on  (3.74)  on  the  left  by  P, 
using  (3.52)  on  the  results  to  get  Pv  =  Pg,  and  then  substituting  g  for  Pg  from  (3.73)  gives 


The  constraint  (3.70)  and 


Pu  =  g. 

Lu  =  LP,ig 


(3.75) 


=  Pug 

=  (Lrr,t)dPg 
=  [Lrfd)dg 
=  Prf , 
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which  follows  from  (3.74),  (2.23),  (2.10),  (3.60),  (3.73),  and  (3.72),  respectively,  show  that  u  is  a  solution 
to  (2.1). 

The  formalism  of  Fishman  and  McCoy  (Ref.  [13,  pp.  47-48])  is  a  subformalism  of  this  report  corre¬ 
sponding  to  (3.70),  ( Lred)d  as  the  reduced  operator  choice,  the  constraint  (3.55), 

Pr  =  P  (3.76) 

for  a  special  class  of  P  and  a  specific  choice  of  L0.  For  this  case,  (3.72)  reduces  to 

{l0  -  P(SL)[L0  +  (I  -  P)8L]-1(5L)}Pu  =  f  (3.77) 

on  use  of  (3.75),  (3.70),  (3.63),  (3.50),  and  (3.76).  Fishman  and  McCoy  use  the  notation  Pv  =  <  v  >  for 
generic  (response  or  stimulus)  v  G  V{L)  otIZ(L),  where  <  v  >  represents  the  macroscale  component  of  v. 
Implicit  use  of  the  admissibility  test  Pf  =  /,  from  (3.70)  and  (3.76),  for  any  given  stimulus  f  is  made  by 
Fishman  and  McCoy  Ref.  [13,  pp.  47-48];  /  is  a  “forcing  with  variations  restricted  to  the  macroscale.”  This 
is  also  true  in  the  formalism  of  Steinberg  and  McCoy  (Ref.  [16,  pp.  1 135-1 136]),  as  seen  by  <  tio  >  =  «o 
in  (19)  and  (22)  of  Ref.  [16,  p.  1135]  when  making  the  notational  association  /  — »  uq.  If  P  is  restricted  so 
as  to  satisfy  the  properties 


P  <A>  =  <  .4  >  P  (3.78) 

PAP  =  <  A  >  P  (3.79) 

for  generic  linear  operator  .4,  and  if  Lq  is  taken  to  be 

Lq  =<  L  >.  (3.80) 


then  PL0  —  L()P,  corresponding  to  (3.48)  with  (3.76),  is  immediately  seen  to  be  true  by  (3.78).  P(8L)P  =  0, 
corresponding  to  (3.55)  with  (3.76),  is  also  true  and  can  be  seen  by 


P{8L)P  =  P{L-  <  L  >)P 

=  PLP  -  P  <L>  P 
=  <  L  >  P-  <  L  >  P2 
=  <  L  >  P-  <  L  >  P 
-  0, 


when  using  (3.49)  and  (3.80),  (3.78),  and  (3.79),  and  P2  =  P,  respectively.  The  constraint  (3.78)  on  P 
corresponds  to  the  constraint  (21)  of  Steinberg  and  McCoy  (Ref.  [16,  p.  1 135]),  but  it  does  not  seem  to  be 
explicitly  acknowledged  by  Fishman  and  McCoy.  The  relation  (3.77)  takes  the  form 

{<!>-<  {SL)[<  L  >  +  (I  —  P)(SL)]-1{8L)  >}<  u  >  =  f  (3.81) 

when  using  (3.80),  (3.79),  and  then  Pit  — ¥  <  u  >.  Equation  (3.81)  is  exactly  the  same  as  the  combination 
of  (5)  and  (6)  of  Fishman  and  McCoy  (Ref.  [13,  p.  48])  if  the  notational  adjustment  SL  — >  L'  is  made  in 
(3.81). 


4.  BISCALE  CONJUGATE  APPROXIMATION 


Consider  a  class  of  methods  under  the  umbrella  of  this  formalism  for  which  the  <  -embedding  for  the  conju¬ 
gate  approximation  is  based  on  the  premise  that  there  are  two  relevant  time  scales,  one  governing  the  “fast” 
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part  of  the  overall  response  and  one  governing  the  “slow”  part  of  the  response.  Specializing  the  derivative- 
expansion  version  (Ref.  [8,  pp.  230-232])  of  the  multiple  scales  perturbation  method  to  a  biscale  expansion 
in  time  prescribes  that 


0  d  d 

d t  dto  +  dt\ 


(4.1) 


be  substituted  for  the  time  derivatives  in  L  to  obtain  T(t ),  where  1 0  is  the  time  variable  associated  with  the 
fast  time  scale  and  t\  is  the  time  variable  associated  with  the  slow  time  scale.  (The  choice  of  two  scales  is 
based  not  only  on  its  simplicity  and  common  occurence,  but  also  that  one  must  be  more  cautious  in  the  use 
of  the  multiple  scales  approach  for  three  or  more  scales  [58].)  This  embedding  leads  to  finite  expansions  for 
T (e)  in  e  for  many  common  V s,  which  are  local  in  time. 


4.1  Frequency  Window  Variant 

A  variant  of  this  temporal  biscale  approach  is  to  use  the  Fourier  representation  of  the  ^-dependence  in  the 
perturbation,  so  that  to -dependence  —>  ^-dependence.  As  to  is  associated  with  the  fast  time  scale,  there 
should  be  a  “cutoff”  value  of  u;0  below  which  the  ti  -dependence  dominates.  For  those  cases  for  which 
L(d/dt)  satisfies 

^W0+eW,i  =  L^  +  ‘W^-  <42) 

where  Jrio  denotes  the  Fourier  transform  of  the  to  -dependence  into  -dependence,  the  continuation  is  taken 
to  be 

T(e<uJo)  =  L(iuJo  +  e  77— ).  .  (4.3) 

at  i 

Substituting  T(e.  lUq)  of  (4.3)  into  (4.2)  leads  to 


(4.4) 

(4.5) 

after  some  rearranging.  This  indicates  that 

£_1  =  3ct0~1T{e,^o)~'i  FioU  =  0-4 1  (4-6) 

can  be  used  in  place  of  (2.33)  as  the  conjugate  approximation  of  L~l .  For  those  cases  satisfying  (4.2), 

To"1  =  T~l(f.^  0)|f_M) 

=  T(wo)"1.  (4.7) 

so  that  To”1/  can  be  interpreted  as  the  time-harmonic  response  of  L  to  /  at  the  frequency  v ,  where 

u/'o  =  2tt7/.  At  a  given  frequency,  the  perturbational  expansion  about  6  =  0  itself  can  be  interpreted  as  a 

transient  departure,  on  the  t\  (slow)  time  scale,  from  the  time  harmonic  solution  at  v  —  The  only 

inversion  required  to  obtain  the  Q j  operators  in  (2.61)  is  T0” 1 ,  so  that  obtaining  the  Qj  ’s  as  a  function  of 
reduces  to  time-harmonic  re-analysis.  A  practical  computational  approximation  is  to  restrict  the  range  of  v 
values  to  a  given  frequency  window  of  interest  for  the  particular  problem  (2.1).  At  least  one  case  satisfying 
(4.2)  is  that  for  which  (2.1)  is  a  set  of  linear  coupled  ordinary  differential  equations  in  time  with  constant 
coefficients,  L  then  is  a  polynomial  in  the  time  derivative  operator. 
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4.1.1  Linear  Coupled  O.D.E.  Systems  with  Constant  Coefficients 

The  ordinary  differential  equation  setting  is  generic  but  at  the  same  time  well-rooted  in  concrete  applications. 
This  setting  encompasses  many  useful  classes  of  models,  with  examples  arising  from  finite  element  methods 
(FEMs)  and  large  linear  electronic  circuit  response.  Additionally,  this  setting  is  sufficiently  specific  in  scope 
as  to  allow  the  explicit  construction  of  the  conjugate  approximation  results  in  terms  of  generic  coefficient 
matrices  for  the  differential  equations.  A  summary  of  these  results  is  given  by  the  following  theorem,  whose 
proof  is  deferred  to  Appendix  E. 

Theorem  9  Let  the  system  operator  L  for  the  Jth  order  set  of  ordinary  differential  equations  (2.1)  be 
denoted  by 

r  4-  r  d  ' 

L  =  Tj-  w» 

3= 0  Ul 

where  the  Ct  coefficients  are  constant,  square  matrices  each  of  the  same  size.  Applying  the  biscale  pertur¬ 
bation  of  (4.3)  to  (4.8)  leads  to  (2.34)  for 


Tj  =  TjM—j 

dtiJ 


TkUo)  =  E 


(w0)m-AX, 


where  the  binomial  coefficients  are 


\  k  J  k\(m  -  k)l 

The  conjugate  approximation  continuation  (2.50)  reduces  to 


(4.10) 


(4.11) 


T(e,u  o)"1  *  V 7, 

j-  0 


where  the  component  matrices  j  are  recursively  given  by 


(4.12) 


=  'Mir1  -  E  k. 


(4.13) 


the  w'o  dependence  being  implied.  The  HjC s  are  defined  by  (2.38). 
The  conjugate  approximation  to  L_1  given  by  (4.6)  reduces  to 


0)^77 


Note  that  (4.9)  and  (4.10)  produce 


Tj  = 


To  =  E  (^0 )mCm 


(4.14) 


(4.15) 


=  To  (4.16) 

as  important  special  cases.  The  bulk  of  the  work  in  obtaining  L~l  is  in  computing  the  time  harmonic 
response  Tq1  for  values  of  ~'o  corresponding  to  the  chosen  frequency  window. 
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4.2  Alternative  ODE-Model  Methods  for  Reduced  Basis 


In  the  case  of  reduced  basis  model  reduction,  an  alternative  variant  of  the  ordinary  differential  equation 
biscale  method  of  Section  4.1.1  serves  as  a  precursor  to  and  a  generalization  of  the  force  derivative  method. 
For  (7  —  L1  L)L~l  —  L~l  -  L1  in  (2.4),  an  alternative  form  for  (2.32)  is 

L -1  *  A(f)|f  =  0— >1  (4-17) 


for 

A  (()  =  LI*  +  \L-l-LI*](e).  (4.18) 

In  this  case,  for  which  the  e -continuation  of  the  conjugate  approximation  is  extended  to  include  the  entire 
residual  term  L~x  -  L1* ,  one  gets 


m 


d)  ^rf  + 


L{dTo+eot1 


dto  dh 


(4.19) 


when  using  reduced  basis  version  (3.14)  of  L1* ,  (4.4),  and  (4.5).  If  <J>,  and  <£»,/  are  independent  of  6  and  u;0, 
then  Lemma  5  of  Appendix  E  can  be  used  with  r  — >•  Qj  — >■  'i/'j,  and  Tj  — >  &r^fj(uJo)$rj  from  (2.34), 
(4.9),  and  T (f)  $r*T (t.  u;0)$rf.  This  results  in 


[QjT{e.>jQ)Qd]-' 


M 


(4.20) 


where  the  component  matrices  t!'?  are  recursively  given  by 


E  Hjk&SToQjr'wSfkQd)  i\j-h 


.A-l 


(4.21) 


the  ~'o  dependence  is  implied.  Substituting  (4.12)  and  (4.20)  into  (4.19)  and  the  results  into  (4.17)  leads  to 

L~l  «  <M$,-t7^r1$rt  + 


M 


dJ 


- <M^o )rt 


to 


(4.22) 


as  a  reduced  basis  hybrid  method  for  sets  of  coupled  ordinary  differential  equations,  where  the  j  ’s  are 
recursively  given  by  (4.13). 


4.2.1  Tunable  Force  Derivative  and  Generalized  Lanczos  Methods 


One  can  generate  a  systematic  approximation  to  the  reduced  basis  hybrid  method  just  developed  by  expand¬ 
ing 


[®j(w’o)  -  $rflV(->o)*rt]  =  E 

A*— 0 


(-■o  -  -•“)*  d[y,  -  *d*'j*S] 


k\ 


<9u 
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about  the  value  lox  as  a  Taylor  series,  assuming  an  analytic  dependence  on  jjq.  If  only  the  first  term  in  the 
expansion  is  kept,  one  gets 


M 


if 


./'=0 


dt  iJ 


*0 


M 


^0_1  E  [*/(<■»*)  -  7 


J= 0 


*0 


i\/ 


M  y 

E  ^(^)^,t]-T 

.7=0 


on  the  right-hand  side  of  (4.22),  so  that 


L~l 


M 


*d(*r'L$d)-l*s  +  E  [$,■(->*)  -  <mE^)<*v+] 


3=  0 


(4.23) 


results  as  the  zeroth-order  approximation  to  (4.22).  It  is  proven  later  that  the  conventional  force  derivative 
method  results  from  the  0  subcase  of  (4.23),  so  that  (4.23)  represents  a  tunable  generalization  of  the 
force  derivative  method. 


Taking  the  above  Taylor  expansion  in  (4.19)  leads  to 

A(0  =  + 

{Tie.^y1  -  (4  24) 

when  keeping  the  first  expansion-term  only  and  canceling  the  Fourier  transform  with  its  inverse.  The  gen¬ 
eralized  Lanczos  reduced  basis  method  is  based  on  the  idea  of  finding  <fv  such  that 

A(e)  =  *d($riL$d)-l$r'  +  0(ek)  (4.25) 

for  (4.24)  for  some  chosen  k  >  0.  For  such  a  choice  of  <f>c/,  only  the  <E>,y  (<T>,. — 1  term  needs  to  be 
kept  in  (4.23).  Such  a  can  be  constructed  if  one  assumes  that  the  admissible  /’ s  to  be  considered  for 
(2.1)  are  to  be  taken  exclusively  from  the  set  V,  where  /  £  V  implies 

f  =  Fg(t)  (4.26) 

for  some  g.  The  matrix  F  common  to  all  of  the  /’ s  is  time-independent,  and  g  is  a  time-dependent  vector.  In 
mechanics  for  example,  (4.26)  represents  a  superposition  of  static  loads  (columns  of  F)  using  time-varying 
coefficients  (components  of  g).  Choose  the  block  columns  of  as  ^  j(^)F  in  the  left-to-right  sequence 
j  —  0. 1 . k  -  1,  so  that 

=  (  §o(~’*)F  $1  (vj“)F  •••  Vk-i(S)F  )•  (4.27) 

where  the  tyj9 s  are  given  by  (4.13).  The  number  of  columns  of  F  (and  correspondingly,  components  of  g) 
should  be  small  compared  to  the  size  of  (4.8),  so  that  the  number  of  columns  of  4>,/,  is  small  compared  to 
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the  number  of  rows  of  <3 Taking  M  -4  k  -  1  and  lJo  — ^  u>*  in  (4.12),  applying  the  result  to  (4.26),  and 
then  substituting  <3^  for  the  right-hand  side  of  (4.27)  leads  to 


T(e.^r\f 


k-l 


( 


J= o  ari 

\ 


f  ) 


t  o 


-llL 


c>U 


,-k-l  Vk-l9 

s//-1 

Ao 


+  0(« 


=  $,/>(<) +  CV) 


for  /?(e)  defined  with  block  rows  f J dJ (j / Oty1  in  the  top-to-bottom  sequence  j  =  0. 1. . . .,  k  —  1.  Using  this 
in 


=  ^[$rtT(f.u;*)#(r1^tr(€,u;>-)r(€,^)-1/+  0(tk) 

=  $rf[^rtr(f,u;-)^]-1*l.tr(€.u>,,‘)^(€)  +  o(ek) 

=  $rffc(0  +  C?(efr) 

=  T(e.^)-1/  +  cV) 

shows  that  (4.24)  reduces  to  (4.25)  for  the  choice  of  (4.27)  for  4>£;.  For  first  or  second  order  (J  =  1  or 
2)  cases,  this  resembles  the  method  of  Haggblad  and  Eriksson  [29].  The  span  of  the  columns  of  for 
4>,/  given  by  (4.27),  will  later  be  shown  to  reduce  to  a  conventional  Krylov  subspace  under  special  circum¬ 
stances.  In  such  cases,  this  becomes  the  conventional  Lanczos  method,  with  proper  normalization  of  4'  /  ’s 
columns. 


4.3  Conventional  Force  Derivative  Lanczos  Submethods 


For  the  linear  dynamic  mechanics  class  of  finite-element  models  for  which  the  force-derivative  and  Lanczos 
methods  were  originally  developed,  the  ordinary  differential  equations  are  second  order  with  real,  constant, 
symmetric  matrix  coefficients.  For  the  symmetric  case,  <E><y  =  $,.  =  <f>  is  reasonable.  The  second-order 
representation  (3.2)  of  (2.1)  is  given  in  terms  of  (4.8)  for ./  =  2  as 

Co  =  I<  (4.28) 

Ci  =  C  (4.29) 

£,  =  M.  (4.30) 


where  u(t)  and  f(t)  of  (2.1)  are  interpreted  as  the  displacement  response  and  applied  loads,  respectively. 
For  an  equivalent  first-order  (or  state)  representation  [38,  p.  12],  L  is  given  by  (4.8)  for  ,7  =  1  and 


Ci 


C-o 


(4.31) 

(4.32) 
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where  the  j£q  and  C]  matrices  are  constant,  real,  and  symmetric  since  this  is  true  for  M,  C,  and  K ,  and 
where 


/ 


u  — y 


are  used  in  (2.1).  The  development  that  follows  will  use  the  second-order  representation. 


(4.33) 

(4.34) 


Equation  (4.10)  with  i*u'Q  — t  oJ  *  reduces  to 

To  =  — u*“M  -j-  4“  Tv 

7)  =  2  +  C 

f2  =  M 


(4.35) 

(4.36) 

(4.37) 


when  using  (4.28)  through  (4.30).  The  reduced  basis  versions  of  these  coefficients  become 

$rf0$  =  -uj*2 1  +  iu>~l3  +  T2 
$rfi$  =  2  iu>*I  +  l3 
<f>TT2<f>  =  I 


(4.38) 

(4.39) 

(4.40) 


when  using  (3.20),  (3.15),  and  (3.21).  The  ihj  component  matrices  are  recursively  given  by  (4.21),  which 
reduces  to 


($Tf0 &)-1 

(4.41) 

(4.42) 

(4.43) 

for  ???  >  2.  The  j  component  matrices  are  recursively  given  by  (4.13),  which  reduces  to 


'J'o 

—  j-1 

—  10 

(4.44) 

=  -T-^T-1 

(4.45) 

'I'm 

=  -To1  +f2*n.-2] 

(4.46) 

for  m  >  2. 


4.3.1  Force  Derivative  as  a  Suhmethod 

The  force-derivative  method  corresponds  to  the  limiting  subcase  of  cj*  — >  0,  so  that  in  the  second-order- 
representation  (4.35)  through  (4.37)  reduce  to 


To  =  I< 

(4.47) 

T \  =  C 

(4.48) 

f2  =  M. 

(4.49) 
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The  reduced  basis  versions  of  these  coefficients  (4.38)  through  (4.40)  become 

4>Tf0$  =  T2 
$tTi<J>  =  3 
<S>Tf2$  =  I. 


The  i'j  component  matrices  of  (4.41)  through  (4.43)  reduce  to 


_  y  —  2 

—  2 


3T~ 2 

['^  t’m—  1 


+  vm-2 


(4.50) 

(4.51) 

(4.52) 


(4.53) 

(4.54) 

(4.55) 


for  m  >  2,  when  using  (4.50)  through  (4.52).  The  'I’j  of  (4.44)  through  (4.46)  reduce  to 


do 

=  A’-1 

(4.56) 

di 

=  —K~XCK~X 

(4.57) 

(4.58) 

for  rn  >  2,  when  using  (4.47)  through  (4.49).  Taking  tv’*  — >■  0  is  somewhat  of  a  contradiction  in  the  context 
of  t0  as  the  fast  variable,  for  which  w’0  (and  hence  w*)  is  associated  with  “higher”  frequencies.  However, 
only  the  / 1  (slow)  time  scale  remains  (take  w’*  —*■  0  and  e  — Y  1  in  (4.3)),  suggesting  / 1  — >  I . 


Using  the  first  three  terms  of  (4.53)  through  (4.55),  the  first  three  terms  of  (4.56)  through  (4.58), 
Y-2  — >  0“2,  3  — >■  A,  w*  — >■  0,  and  ty  — >  t  for  the  slow  time  all  in  (4.23)  reproduces  (38)  of  Ref.  [38, 
p.  29].  (The  Y-2  and  3  substitutions  account  for  notational  discrepencies  between  Ref.  [38]  and  this  re¬ 
port.)  For  the  case  of  an  arbitrary  number  of  terms,  taking  t !  — t  t  and  uj*  — >  0  in  (4.23),  substituting  this 
result  for  L-1  into  A-1/  (from  (2.1)),  and  then  using  the  notational  associations 


6  < — ^  $ 

(4.59) 

q  «— >■  (4>tA$)-1$t/ 

(4.60) 

f)j 

^  i — >  - r  f  for  j  >  0 

dt \:r 

(4.61) 

i.j  < — »  t'j  for  j  >0 

(4.62) 

i.j  < — >  d ,  for  j  >  0 

(4.63) 

reproduces  the  conventional  force  derivative  results  (39)  of  Ref.  [38,  p.  30].  Proof  of  this  essentially  hinges 
on  proving  the  correspondence  between  the  ,  and  'f ■  coefficients  of  (4.53)  through  (4.58)  of  this  report 
with  the  .4ij  and  By.,  coefficients  of  Ref.  [38,  p.  30]  and  [59,  p.  716].  To  show  this,  the  equations  given 
by  both  ReT  [38,  p.  30],  and  by  (A8),  (A9),  (A16),  and  (A17)  of  Ref.  [59,  p.  716],  for  the  .4,.,  and  By., 
coefficients  are  reproduced  here  as 


-4ij  =  -H-2A.41.,_1  -  ft  21,h  for  j  >  1 
■h.j  =  -4i.,_i  for  j  >  1 

.i,.o  =  ft-2 
-4-2.0  =  o 


(4.64) 

(4.65) 

(4.66) 

(4.67) 
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and 


fiij  =  for  j  >  1  (4.68) 

Bi.j  =  forj  >  1  (4.69) 

Buo  =  A'"1  (4.70) 

B2.o  =  0.  (4.71) 

respectively.  (The  notation  ( 1  2  and  A  of  Ref.  [59]  corresponds  to  T-2  and  3  of  this  report,  respectively. 

The  “hats”  on  Q~2  and  A  of  Ref.  [38]  have  been  “dropped”  in  favor  of  the  notation  of  Ref.  [59].)  It  is  readily 

verified  that  T-2  and  A  — »  3  in  (4.66),  (4.67),  and  the  j  =  1  case  of  (4.64)  lead  to 

Aj.!  =  -Y~2/3T-2,  (4.72) 

and  that  (4.70),  (4.71),  and  the  j  —  1  case  of  (4.68)  lead  to 

fii.!  =  -K~lCK~l .  (4.73) 


Equations  (4.53),  (4.54),  and  (4.62)  are  in  agreement  with  ft-2  >  T-2  in  (4.66)  and  (4.72),  and  equations 
(4.56),  (4.57),  and  (4.63)  are  in  agreement  with  (4.70)  and  (4.73).  Equations  (4.62)  and  (4.63)  are  verified 
for  j  =  0.1.  Taking  j  — >  j'  —  1  and  dropping  the  prime,  equation  (4.65)  leads  to 

,42j-i  =  .4ij_2  forj  >  2. 

and  (4.69)  leads  to 

S2lj-i  ~  forj  >  2. 

Substituting  these  into  (4.64)  and  (4.68)  leads  to 

-4ij  =  -T -  T-2.4l!7_2  forj  >  2  (4.74) 

and 

Bhj  =  -K~lCB\.j-i  -  K_1  MBij-2  for  j  >  2,  (4.75) 

respectively,  when  using  Q~2  -r  T-2  and  A  ->  0.  Substituting  (4.62)  into  (4.74)  and  (4.63)  into  (4.75) 
leads  to  (4.55)  and  (4.58),  respectively.  Equations  (4.62)  and  (4.63)  are  verified  for  j  >  2  as  well.  Taking 
ti  -a  t  and  ^  0  into  (4.23),  substituting  this  result  for  L~l  into  T-1/,  and  then  using  (4.53)  through 
(4.58)  in  the  results  is  collectively  equivalent  to  the  conventional  force  derivative  results  of  Ref.  [38,  p.  30] 
for  second-order  formulated  systems. 

43.2  Lanczos  as  a  Submethod 

In  the  limiting  subcase  ^  — y  0,  (4.13)  reduces  to  (4.56)  through  (4.58)  for  the  T^’s  defining  $  in  (4.27).  If 
one  takes  C  —  0  in  (4.56)  through  (4.58),  they  further  reduce  to 

=  A'-1  (4.76) 

(i>2/— i  =  0  for  j  >  1  (4.77) 

=  (— A'-1J\/)J  A’-1  for  j  >  1,  (4.78) 

so  that  (4.27)  reduces  to 

4>=(a-1F  (-A-1M)A'-,A  •••  (-K-WIYK-'f) 


(4.79) 
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after  eliminating  the  zero  columns  arising  from  the  T'j ’s  of  odd  j,  where  (  =  (k  -  l)/2  for  even  k  —  1  and 
C  =  ( k  -  2)/2  for  odd  £  -  1.  Except  for  the  (irrelevant)  alternating  sign  with  each  power  of  -K~1M,  the 
columns  of  3>  in  (4.79)  are  seen  to  form  a  block  version  of  the  Krylov  sequence  (Ref.  [24,  p.  566])  tradition¬ 
ally  associated  with  L  specialized  to  J  =  2  in  (4.8).  The  initial  block  A  _1F  of  (4.79),  whose  columns  are 
the  static  solutions  to  the  loads  of  the  columns  of  F,  is  the  usual  recommended  initializing  choice  (Ref.  [24, 
p.  568])  for  the  Krylov  sequence  in  reduced  basis  applications.  (The  case  for  which  F  has  just  one  column 
and  (j  just  one  component  results  in  the  standard,  “nonblock”  version  of  the  Krylov  sequence.)  The  con¬ 
ventional  Lanczos  method  uses  a  reduced  basis  consisting  of  the  above  Krylov  sequence,  which  has  been 
orthonormalized  with  respect  to  the  mass  matrix. 

One  might  suspect  that  the  C  =  0  assumption  may  hurt  the  Lanczos  method’s  ability  to  handle  general 
damping  cases,  and  some  evidence  of  this  is  found  in  the  conclusions  of  Ref.  [36].  In  particular,  the  force- 
derivative  method,  which  uses  a  finite  set  of  the  nonzero-C  4'/ s  given  by  (4.56)  through  (4.58),  efficiently 
and  accurately  handles  the  nonproportional  damping  case  of  Ref.  [36],  in  contrast  to  the  Lanczos  method. 
Nevertheless,  one  possible  advantage  to  the  usual  Krylov  sequence  (from  the  C  =  0  assumption)  in  con¬ 
junction  with  the  orthonormalization  process  is  that  together  they  produce  a  tridiagonal  reduced  problem 
(see  Nour-Omid  and  Clough  (Ref.  [24,  pp.  567,  569])  or  Golub  and  Van  Loan  (Ref.  [60,  p.  477])). 
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Appendix  A 


A.  PROOF  OF  THEOREM  2 

The  following  lemmas  are  preliminary  to  proving  theorem  2. 

Lemma  1  Given  operators  Pd,  Pr,  and  T,  the  operator  Pd(PrT Pj)1  Pr,  when  it  exists,  is  an  outer  gener¬ 
alized  inverse  of  T. 

Proof:  Substituting  Pd(PrT  Pd)1  P,-  for  A1  in  A1  A  A1  leads  to 

PdiPrTPdyPrTPdiPrTPdyPr  =  Pd(PrT  Pd)1  (PrT  Pd){PrT  Pd)1  Pr 

=  PdiPrTPdYPr , 


proving  the  lemma. 

Lemma  2  Given  a  fixed  linear  operator  T  and  two  fixed  linear,  idempotent  Pd  and  Pr ,  PdT1  Pr  is  unique 
for  all  Tl  such  that 

PdT !T  =  Pd  (A.l) 

TT1?,.  =  P,.  (A.2) 

are  simultaneously  true. 

Proof:  Suppose  that  there  are  two  outer  generalized  inverses  of  T,  label  them  T/  and  T2 7,  each  of  which 
satisfied  (A.l)  and  (A.2)  simultaneously.  They  would  obviously  satisfy 

PdTiLT  =  PdT-jT  (A.3) 

TTdp.  =  TT/Pr,  (A  .4) 


and  because  they  are  outer  generalized  inverses  of  T,  they  would  also  satisfy 


and  hence 


r/TT/  =  rd 
TdTPj  =  ToJ, 


PdTdTTdp.  =  PdTdP 
PdT/TTdPr  =  PjTl'Pr- 


Substituting  (A.4)  into  (A.5),  and  (A.3)  into  (A.6),  leads  to 

PdTdTTyp  =  FVr/P, 
PlTdjP'Pr  =  P.IP'P.. 


(A.5) 

(A.6) 


respectively.  These  final  two  equations  clearly  show  that  PjTi1  Pr  =  PdT}1  P,  ,  and  the  lemma’s  proof  is 
complete. 
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Lemma  3  Given  two  linear  idempotent  operators  P  and  II, 

PU  =  U 
72(11)  C  72(F) 
PU  =  P 
7 2(F)  C  72(11)  . 


(A.7) 

(A.8) 


Proof:  If  /  £  72.(11 )  then  t  =  Ut  =  PUi  —  Pt  so  that  t  £  72(F),  proving  (A.7).  The  relation  FIT  =  P  implies 
P(I  -  II)  =  0,  which  implies  (I  -  F)  (I  -  II)  =  (I  -  II),  which  in  turn  implies  72(7  —  II)  C  72(7  -  P)  by 
(A.7).  This,  in  conjunction  with  72(7  —  II)  =  A"(II)  and  72(7  —  P)  =  A  '(F),  leads  to  A  (II)  C  A "(F).  Us¬ 
ing  the  fact  that  .4 c  C  Bc  impliesF  C  .4  forgeneric  .4  and  B,  A^Il)6  =  72(II)-{0},andA'(F)<r  =7 2(F)- 
{0}  in  the  results  leads  to  72(F)  -  {0}  C  72.(11)  -  {0},  where  the  superscript  C  denotes  the  (set)  comple¬ 
ment.  The  union  of  both  sides  of  this  result  with  {0},  along  with  0  £  72(11)  and  0  £  7 2(F),  leads  to  (A.8). 

The  proof  of  Theorem  2  itself  is  now  given  as  follows. 

•  To  show  that  L1'  is  an  outer  generalized  inverse  of  L,  use  (2.10)  to  substitute  PrLPd  for  Le/j  in 
(2.1 1),  and  then  apply  Lemma  1,  with  T  — >  L,  to  the  result. 

•  To  prove  (2.12),  let  v  £  A'(  PrL).  One  then  has  (7  —  L1  L)v  —  (I  —  F/ L f  j jJ PrL)v  =  r  from  (2.1 1) 
for  idempotent  (I  -  L1* L),  so  that  v  £  TZ(I  —  Lr  L),  proving  (2.12). 

•  The  relation  L1"  LPj  =  PriLfj  f1  PrLP,j  =  P(iLej  j1  Lef  j  =  Fj  follows  from  (2.11),  (2.10), and  (2.13), 
proving  that  (2.16)  follows  from  (2.13).  The  converse,  that  (2.13)  follows  from  (2.16),  is  also  clearly 
true. 

•  To  prove  (2.15)  follows  from  (2.13),  one  first  notes  that  application  of  (A.7)  of  Lemma  3 ,  with  F  — >•  Pj 
and  n  ->  LrL,  to  PdLr L  =  Pd2LFjSlPrL  =  PdL.j/p.L  =  LrL  leads  to  1Z(Lr L)  C  7 Z(Pd). 
However,  application  of  (A.7)  of  Lemma  3,  with  F  — >  L1* L  and  II  — >  Pd,  to  (2.16)  leads  to 
R{Pd)  F  'R[Lr' L).  Together  these  results  prove  F(L^L)  =  lZ(Pd),  which  combines  with 
■R{LIXL)  =  A •'(/  -  Lr*L)  to  prove  (2.15). 

•  To  prove  (2.17)  followsfrom  (2.13),  one  starts  with  [R(Pd)  r\M{Pr  L)]  C  [; R(Pd )  r\R(I  —  L1  L)]  = 

| A"( /  —  L1* L)  fl  R(I  —  L1*  L)]  =  {0},  which  follows  from  (2.12)  and  (2.15).  Substituting  this  into 
A '{Lf  j j)  =  {[R(Pd)  fl  X{PrL)]  U  .V'(F/)},  which  follows  from  (2.10),  leads  to  (2.17)  upon  using 
0  £  A'(F/). 

•  The  relation  PrLLr  =  P,-LPdLf  } / P,  =  L.jfL^/P,.  =  F,  followsfrom  (2.1 1),  (2.10),  and  (2.14), 
proving  that  (2.19)  follows  from  (2.14).  The  converse,  that  (2.14)  follows  from  (2.19),  is  also  clearly 
true. 

•  To  prove  (2.18)  follows  from  (2.14),  one  first  notes  that  application  of  (A.8)  of  Lemma  3,  with 
F  LLr  and  II  -4  Fr,  to  LlX Pr  =  LPdL.j/P2  =  LPdLtSSlPr  =  LLr  leads  to 
R(LL1*)  C  72(F  ).  However,  application  of  (A.8)  of  Lemma  3,  with  F  — »•  Pr  and  II  — >  LL I~,  to 

(2.19)  leads  to  R(Pr)  C  R(LlJ~).  Together  these  results  prove  R(LLr)  =  R  (F,),  which  combines 
with  R(LLI>')  =  A  '{I  —  LL1")  to  prove  (2.18). 

•  To  prove  (2.20)  follows  from  (2.14),  firstnote  that  F,L(/ —  =  (Pr  —  PrLLl’‘)L  =  (. Pr  —  P,)L  = 

0  follows  from  (2.19),  so  that  72(7  -  L1" L)  C  A'(FrL).  This  result  combines  with  (2.12)  to  prove 

(2.20) . 
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•  As  T  Lejf  in  (A.l)  and  (A.2)  of  Lemma  2  corresponds  to  (2.13)  and  (2.14),  respectively,  then 
(2.11)  and  Lemma  2  together  show  that  L 7*  is  unique  when  (2.13)  and  (2.14)  are  simultaneously 
satisfied. 

•  The  relation  Lr  L  =  PdLej j1  P,  L  =  Pd.Ltj  j1  PrLPc{  =  PdLej j1  Lejj  =  P{ j  follows  from  (2.11), 
(2.21),  (2.10),  and  (2.13),  proving  that  (2.22)  follows  from  (2.13)  and  (2.21). 

•  The  relations  (2.21)  and  (2.22)  imply  0  =  PrL(I  -  Pd)  =  PrL(I  -  LrL)  =  (Pr  -  PrLL^)L,  so 
that  7 Z(L)  C  Ar{Pr  —  PrLLr).  This  proves  (2.16),  from  which  (2.13)  follows. 

•  The  relation  LL1"  =  LPjLef j1  Pr  =  PrLPdLeffIPr  =  Lej jLej f1  Pr  ~  Pr  follows  from  (2.11), 
(2.23),  (2.10),  and  (2.14),  proving  that  (2.24)  follows  from  (2.14)  and  (2.23). 

•  The  relations  (2.23)  and  (2.24)  imply  0  =  (/  -  Pr)LPd  =  (/  -  LLr)LPd  =  L{Pd  -  LrLPd),  so 
that  lZ(Pd  -  L!'  LPd)  C  N{L)  =  {0}.  This  proves  (2.19),  from  which  (2.14)  follows. 


The  proof  of  Theorem  2  is  complete. 
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(B.2) 

(B.3) 


(B.4) 

(B.5) 

(B.6) 
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+  E  <J  [(e  -  //  v.,r/?, 

j=M+ 1  L\a-=o  / 

=  X>'  [  (E  -  //.v,r  ^ 

.7=0  L  \k=0  / 

+  E  -HxjTRj 

j=M+ 1  L  \k=0  / 

j=0  LVa— 0  / 

Max(A'..v+J)  /  j  \ 

+  E  E  HjkHALj.k.TkQRj.k  -  - 

j=A/+l  \A-=0  / 

The  step  from  the  second  to  the  third  line  uses  HM.j-k  =  1  for  j  <  AI.  The  last  line  results  from 

j  j  >  Max(A\  M  +  J)  and  1  f  a-  =  0  and  1 

\0<A  <J  j  \  H\j  =  0  J- 

which  follows  from  { j  >  Max(A\  il/  +  •/)}—>■  {j  >  A/  +  J  and  j  >  A  }, 

so  that  H\ j  =  0  and  j  —  k  >  i\/  +  J  —  A*;  but  {0  A  k  <  J}  — >  j:'!/  +  J  —  k  >  3/},  so  that  j  —  A*  >  j\/ 
and  hence  Fl^p.^k  —  0. 

Finally,  the  previous  result  establishes  the  equivalence  of 

E  f  ( E  -  hNj rRj  =  o 

.7=0  L  \fc= 0  / 

and  (2.41).  Solving  this  for  T0QR  j  leads  to  (2.36).  The  proof  of  (2.37)  and  (2.42)  from  (B.6)  is  essentially 
the  same  except  that  QLj-^T a  and  I  1‘ ;  replace  Tk£lRj-k  and  r^,,  respectively. 


Appendix  C 


C.  PROOF  OF  THEOREM  5 


One  can  prove  QR  j  =-  QL  (  for  each  j  >  0  by  induction,  where  the  initial  case  for  j  =  0,  QRo  =  0/'n  =  -4, 
follows  directly  from  (2.45)  and  (2.46).  The  strategy  is  to  first  assume  that 

fiV*  =  OPj-k  (C.l) 

for  !</,<  min  (J.  j)  and  then  prove  that  (2.45)  and  (2.46)  subsequently  lead  to  QRj  =  QLj.  Substituting 
(C.l)  into  the  left-hand  side  of  (2.46)  with  m  — >  j  —  k  leads  to 


=  Hoj-kA  -  V  Hj-M&j-ic+vT'A. 


C=  1 

Substituting  this  result  into  the  right-hand  side  of  (2.45)  with  m  —¥  j  leads  to 

J 

VR,  =  H0jA-Y,HjkATkQRj-k 
k=  1 
J 

=  H0jA-Y,HjkATk 


k—1 

J 


e=i 


Hoj-kA-'£Hj-k,caLj_(e+k)TtA 
=  -  Y  HjkHoj-kATkA 


k~\ 


J  J 


+E  E  HJkHj-k.tATknLj-lc+lf)T(A 

A-l f=i 

J 

H0jA  -  J2  HjcHo.j-eAT(A 


(=1 


j  j 


CEw- 

(=1  A-=l 

-c.kATk$lL  j_({+k)TcA 

J 

j 

=  H0JA  V  HjC 

Ho.j-cA  -  Y  Hj-ckATk&j-it+k) 

/  =  1 

k= 1 

J 

J 

=  HojA~YH.it 

(~X 

Hoj-cA  -  Y  //  -c.kATk.QRu_()_k 
A-l 

TiA 


TCA 


H0jA  -  HjcQ^j-tTiA 


(  =  1 

J 


HojA  -  Y  HjcQLj-(T(A 


C= 1 
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when  using  HjkHj-k.i  =  //,v  a-  (proof  follows),  tlL(j_c)-k  =  by  (C.l);  (2.45)  with  m  ->•  j- 

(,  j~(  =  &Lj-t  by  (C.l)  since  1  <  (  <  min (J,j)  in  the  summation  for  C,  and  (2.46)  with  m  —r  j,  re¬ 

spectively.  The  proof  reduces  to  that  of  proving  that  HjkHj-k.c  =  //  //  -  ./,  for  1  <  k  <  min  (•/.,/)  and 
1  <  (  <  min (J.j). 


The  general  property 


for  all  integers  i,  j,  and  k,  leads  to 


Hi±k.j±k  =  Hjj 


Hj-k.(  —  Hjx+k 
=  Hj-i.k- 


(C.2) 


One  can  prove 


HjkHjj+k  =  H.i.c+k 


(C.3) 


as  follows:  one  has  either  Hjj+k  =  0  or  Hjx+k  =  1.  An  assumption  of  HhL+k  =  0  implies  that 
HjkHj.c+k  =  0.  An  assumption  of  HjA+k  =  1  implies  k  +  C  <  j,  which  in  turn  implies  k  <  j  because 
all  values  of  (  are  positive,  so  that  H,k  =  1  and  HjkHj.c+k  =  1-  As  (  and  k  have  exactly  the  same  range  of 
values,  the  proof  of  the  relation 

HjcHj.e+k  =  Hjj+k  (C.4) 

is  essentially  the  same  as  that  of  (C.3).  In  fact,  (C.4)  is  just  a  relabeling  of  (C.3).  Substitute  I from 
(C.2)  into  the  left-hand  side  of  (C.4),  substitute  Hj-k.c  from  (C.2)  into  the  left-hand  side  of  (C.3),  and  equate 
the  results  (the  right-hand  sides  of  each  are  equal)  to  prove  HjkH  -  HjtHj-c.k ■ 


Appendix  D 


D.  PROOF  OF  THEOREMS  6  and  7 


Proof  of  Theorem  6:  The  proof  of  (3.29)  is  given  by 

Pr$r  =  V~lXlaV*r 

=  v~ln  L,#r 

v-\  (  Irani  0  )  f  r]m  \ 

"  V  «  0  ){  \*r]'  ) 

-p-l(  \ 

V  «&-]m  J 

_  V-1  (  [*r]m  \ 

V  [*r],  ) 

=  V~l^r 

=  $r, 


using  (3.35),  (3.32),  y  — >■  a  in  (3.33),  (3.31),  (3.34),  (3.31)  again,  and  (3.32)  again,  respectively.  The 
preliminary  result 


V-'LYlnV  = 


/  ^  mm 

ms  | 

f  *mm 

V  L$m 

Lss  I 

ft 

o  J 

V 


=  p-'n^ev 

follows  from  using  (3.28)  and  y  — >  ;i  in  (3.33),  Lsm  +  Lssf  =  a[Lmm  +  Lms3]  from  (3.36),  (3.38),  y  — >■  a 
in  (3.33),  and  (3.40),  respectively.  Substituting  LPj  —  P~1LPP~1Ua,P  =  P~lLY\  iP  from  (3.27)  and 
(3.37)  into  the  previous  result  leads  to 

LPd  =  P~lnoOP.  (D.l) 
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Matrix  multiplying  (D.l)  from  the  left  by  F,  and  substituting  Pr'P~lY\0  =  P~lU0PP~lU0  =  P_1no2  = 
T7*- 1  no ,  from  (3.35)  into  the  result  leads  to  LPd  =  PrLPd,  which  proves  (2.23).  This  gives  LPd  =  Lf  jj, 
when  using  (2.10),  which  combines  with  (D.l)  to  give 

LfJJ  =  T~1  n,,0P.  (D.2) 

The  Le  f  j1  of  (3.39)  is  proven  to  be  an  outer  generalized  inverse  of  Lf  j j  from  (D.2)  by  Lf_j j1  Lf  j fLej j1  = 
F_1@7no@@7F  =  F-J070@7'P  =  P_107P  =  Z,e//7,  where  07ri,  =  07  is  easily  verified  from  -(  -4  a 
in  (3.33)  and  (3.41).  The  relation  (2.13)  is  proven  by  PdLr:j  j1  Lejj  =  P-1n,307no0P  =  P_1n,307©'P  = 
V-lY[sW0V  =  T~lYUT  =  Pd  when  using  07no  =  07  again;  07©  =  ll0  from  y  -4  0  in  (3.33),  and 
n.3n0  =  II  i-  The  relation  (2.14)  follows  from  Lef =  P_1no007  F  =  P-1II,  IIoF  =  P_1noP  = 
Pr.  Finally,  using  ©7no  =  07  again  in  L1*  =  F(/Z/f//7F,.  =  F_1n,3©7fl.,  F  proves  (3.42). 

Proof  of  Theorem  7:  Matrix  multiplying  (D.  1 )  on  the  right  by  P~Pffhf\m  0 )  ,  substituting  for  0  from 
(3.40),  substituting  for  F(/F-1  =  P-1n  3  from  (3.37),  and  then  substituting  from  (3.45)  gives 


p~l  n„ 


Lred\^d\r 

0 


=  LP~lU3 


=  LQ,i. 


Right  matrix  multiplying  (3.44)  by  ([4>r],„  0)  gives 

r-1  n,,  =  <*>,(  o  )• 

when  using  ([<f>,.]n,T  0)r([4>,.]ni_1  0)  =  n0  for  ^  -4  0  in  (3.33),  and  non0  =  II,.  Substituting  this  for 
V~l  no  into  the  previous  result  and  then  matrix  multiplying  on  the  left  by  leads  to 


L  red[^d]m 

0 


=  (*,•*$,)(  [<J>,]m  1  0 

=  (QMrftQrU-1  Lrfd[*d]m. 

Inverting  (D.3)  and  matrix  multiplying  the  result  on  the  left  by  <3>r/  gives 

=  Qd[Qd]m-1Lrfd-1[&r]n,{&S&r)-'1 

=  p-1  n..^  7,Jm 
=  p~1  n.^  Lred 

_  tt-Io  (  Lpfrf  *  0  V-p^-l  (  I  mm  0  \  [ 


T  nf  0  0 

p-lihe'pp-'  n„  ( 


pp- 

! 

&■]■ 

o 


a  0  j 
(4>,t4>r)~1 


//^p^ri,  (  ^0]"‘  )  (4>rt<j>l.)_1 


Lrsub<f>,.(<bM,.)-', 


(D.3) 


A  Generic  Bilevel  Formalism  for  Unifying  and  Extending  Model  Reduction  Methods 


47 


when  using 


*d[*d]m  1  =  r-'n^  Imf  j 

from  (3.45),  and  then  (3.41),  ^  — y  a  in  (3.33),  (3.42),  and  (3.44),  respectively.  Matrix  multiplying  on  the 
right  by  and  then  substituting  for  Lrrb  from  (3.14)  leads  to  (3.46). 
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E.  PROOF  OF  THEOREM  9 

The  proof  consists  of  the  following  two  lemmas:  Lemma  4  applies  directly  to  the  theorem  with  no 
change,  and  Lemma  5  applies  with  r  -¥  t 1  and  U,  -4  T;. 

Lemma  4  Under  the  hypothesis  of  Theorem  9,  applying  the  biscale  perturbation  of  (4.3)  to  (4.8)  leads  to 
(2.34)  for  (4.9)  and  (4.10). 

Proof:  Substituting  (4.8)  into  (4.3)  gives 


T(e. 


xy(  i  )(-o r* 

k=0 
J 


k  dt s 

i;  ".•••'(  l  W«r‘- 


k=0 

J 


J 


^o)  ~  E  CjjiuQ  +  ed/ 8ti]J 

j= o 
J 

= 

i=o 

J 

-  £*i 

o 

= 

A-=0 
J 

- 

A— 0 
J 

= 

A'=0 

=  X>fcr* 

A— 0 


0* 


a* 


/ 


L>o)J-kCj 


54  u 

dtik 


dt!h 
dk 


dti 


d^_ 

duk 


when  using  a  binomial  expansion  (exact  because  it  is  finite),  the  definition  (2.38)  of  Hj\ 

■  twice,  (4.10),  and 

(4.9),  respectively. 

Lemma  5  If  Q;  satisfies  (2.61)  and 

Tj 

'*>  ft 

II 

(E.l) 

Elf. 

()Tk  J 

f  dk 

~ 

(E.2) 

for  0  <  k  <  oc  for  each  j  of  0  <  j  <  J,  then 

a, 

=  1 
c>r' 

(E.3) 

=  44 

(E.4) 
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for  0  <  k  <  oc  for  each  j  of  0  <  j  <  M,  where  the  ft ;  component  operators  are  recursively  given  by 


=  H0jfo'  -  E  Hjkfo'TkQj-k.. 

A=1 


(E.5) 


Proof:  Equation  (E.l)  for  j  =  0  gives  T0  =  To  and  T0  1  =  T0  1 .  Operating  on  both  sides  of  (E.2)  with 


j  =  0  by  T0  1  leads  to 


f-i  _  f-i 
drk  °  °  0rk' 


(E.6) 


which  shows  that  (E.4)  is  true  for  j  =  0  when  using  fto  =  from  (E.5).  For  j  =  0,  equation  (2.61)  gives 
Q0  =  T0~l ,  which  combines  with  Q0  =  ff 1  and  T0~l  =  Tf 1  to  give  £20  =  ^o,  in  agreement  with  (E.3)  for 
j  =  0.  Equation  (E.3)  is  also  true  for  j  =  0.  To  prove  (E.4)  by  induction,  take  j  — >  m  in  (E.5)  for  m  >  0 
and  operate  from  the  left  by  the  linear  operator  dk /i hk  to  get 


dk 

drk 


Q, 


i)k  J  ^  _ 

—(-(Z^kT^TkQm-k]) 


dr 


k.=i 

J  >)k 


k=i  °t 

j  ,  -  if  ~ 

-[YH^kTo1T,.-Tnm-k] 

k= 1 

J  .  _  Qk 

-tE  Hmk.TflTk.nm-k}—1: 

A— J  C'1 


a, 


dk 

9rk 


when  using  (E.6),  (E.2)  for  j  =  k,  (E.4)  for  all  j  <  m  -  1,  and  j  -4  rn  in  (E.5)  for  rn  >  0  again,  respectively. 
Equation  (E.4)  is  proven  by  induction  if  the  il ,  component  operators  are  given  by  (E.5).  To  prove  (E.3)  by 
induction,  assume  that  (E.3)  is  true  for  all  j  <  m  —  1  (proven  for  j  =  0)  for  some  finite  in  >  0.  Equation 
(2.61)  and  m  >  0  give 


-V  HmkT0-lTkQm-k 


k=l 

J 


■:ym  —  k 


d 


=  — E  Huik-Tu-1  TkQm-k  — 


0rk  "c h"l~k 

dm 


=  ft 


—k 

k=\ 

J  _  Qk  _ 

-E 

k—l 
J 

-E  HmkT0-lfkQm 

A-=l 
J 

-[E  Hnkfo'na 

A-l 

-  <9m 


-At 


()Tm 

dm 


dTm 
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when  using  (E.3)  for  all  j  <  m  -  1,  (E.l),  (E.4)  for  all  j  <  m  —  1,  To  1  =  T0  1,  and  (E.5)  for  m  >  0, 
respectively.  Equation  (E.3)  is  also  proven  by  induction  if  the  Qj  component  operators  are  given  by  (E.5). 


